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SUMMARY 


A  detailed  description  is  given  of  a  program  which  has  been  developed 
to  study  analytically  the  performance  of  a  UTIAS  Implosion-Driven  Shock 
Tube.  Some  results  of  the  computations  are  presented.  These  have  been 
checked  experimentally  using  an  existing  8  inch  diameter  implosion  chamber 
and  a  5/16  inch  diameter  channel.  The  experimentally  measured  shock-wave 
Mach  numbers  were  much  greater  than  anticipated  from  the  computations.  The 
cause  of  the  discrepencies  is  shown  to  be  largely  computational,  and  reco¬ 
mmendations  to  improve  the  calculations  are  proposed.  The  experimentally 
measured  shock  Mach  numbers  decrease  quite  rapidly  in  the  5/16  inch  diameter 
tube,  as  expected.  Nevertheless,  it  was  possible  to  obtain  a  shock  Mach 
number  M,.  -  **0  in  air  using  only  a  600  psi  2Hg  +  O2  gaseous  implosion  driver. 

Much  higher  shock  Mach  numbers  are  expected  with  a  coupled  PBTN  explosive 
driver . 
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1 .  INTRODUCTION 


In  order  to  overcome  the  current  "performance  barrier"  of  Hypervelocity 
Launchers,  Professor  I.  I.  Glass  proposed  in  1959  the  use  of  spherical,  im¬ 
ploding  shock  waves.  Since  then,  the  Institute  for  Aerospace  Studies  has  been 
actively  engaged  in  this  field  1.  An  Implosion-Driven  Hypervelocity  Launcher 
has  been  designed  and  Fig.l  illustrates  schematically  its  principle  of  operation. 
A  model  with  an  8  inch  diameter  hemispherical  chamber  has  been  built  and 
several  experimental  studies  are  going  on.  A  numerical  model  has  also  been 
developed  and  some  previous  investigations  2-3  have  been  concerned  with  the 
computed  performance  and  the  optimization  of  this  launcher. 

The  present  preliminary  study  deals  mainly  with  the  computation  of  the 
performance  of  an  implosion-driven  shock  tube.  The  primary  aim  was  to  study 
analytically  the  limit  of  a  launcher  as  the  mass  of  the  projectile  goes  to 
zero.  Such  a  study  should  enable  us  to  place  an  upper  limit  on  the  projectile 
velocities  that  could  be  attained  with  this  type  of  hypervelocity  launcher. 

It  was  thought  also,  following  the  first  results,  that  an  implosion-driven 
shock  tube  would  be  a  good  facility  to  obtain  extremely  high  shock-wave  Mach 
numbers  at  relatively  high  channel  pressures. 

An  analytical  investigation  is  presented  here  and  is  followed  by  some 
experimental  results  obtained  in  a  5/l6  i&ch  diameter  channel  (barrel)  with  an 
8  inch  diameter  hemispherical  chamber.  This  experimental  work  was  aimed  at 
checking  the  results  of  the  computations  and  to  obtain  an  estimate  of  the 
losses,  as  the  numerical  code  is  based  on  an  inviscid,  adiabatic  flow  field. 

As  it  will  be  seen,  the  agreement  between  experimental  and  computed  results 
was  found  to  be  poor.  The  reasons  for  these  differences  are  shown  to  result 
from  computational  limitations  and  recommendations  for  future  work  on  this 
subject  are  given.  The  actual  measured  large  shock  wave  attenuation  (for 
example,  using  a  600  psi  2Hg  +  Og  driver  and  a  5/16  inch  diameter  channel  the 
shock  Mach  number  decreases  Ms  ~  40  to  Mg  ~  10  in  JO  inches  or  AMg  =  5/ft,  see 
Fig. 21 )  can  be  attributed  to  viscous  effects  primarily,  and  will  be  investigated 
in  greater  detail.  The  5/l6  inch  diameter  barrel  is  not  a  realistic  shock- tube 
channel.  As  noted,  it  was  used  to  check  the  limiting  case  of  a  massless  pro¬ 
jectile.  In  future,  1  inch  diameter  channels  will  also  be  utilized  for  this 
study. 

2.  THEORETICAL  CONSIDERATIONS 


In  this  section,  a  complete  description  of  the  program  that  has  been  used 
in  this  study  and  its  theoretical  background  are  given.  This  code  was  developed 
independently  by  Piacesi  14  and  Sevray  It  has  been  used  by  Sevray  to  analyze 
the  operation  of  the  UTIAS  Hypervelocity  Launcher  and  to  optimize  its  perfor¬ 
mance.  It  was  also  used  with  some  modifications  by  Flagg  3  in  order  to  design 
a  second  generation  launcher,  UI1AS  Implosion-Driven  Hypervelocity  Lhanhhet 
Mark  II.  Although  these  calculations  have  been  very  useful  in  predicting  the 
analytical  performance  of  the  launcher,  the  actual  experimental  performance 
is  only  of  the  order  of  h<y£  of  the  computed  performance,  consequently,  it  will 
be  necessary  to  take  account  of  the  effects  of  radiative,  convective,  ablative 
and  frictional  losses  to  make  the  calculations  more  realistic. 

For  the  present  work,  the  code  has  r.>een  auopted  to  handle  the  shock  tube 


case  essentially  by  adding  a  third  region  (Fig. 2)  in  front  of  the  diaphragm  or 
of  the  projectile,  which  also  permits  one  to  analyze  theoretically  the  influence 
of  counterpressure  on  the  velocity  of  the  projectile. 

2.1  Description  of  the  Code 

This  numerical  code,  written  in  a  Lagrangian  form,  is  based  upon  the 
artificial  viscosity  technique'  as  established  by  von  Neumann  and  Richtmyer  5 
to  handle  shock- wave  problems. 

The  system  (Fig. 2)  is  divided  into  three  regions:  the  explosive  PETN,  the 
gas  mixture  2Hp  +  Op  and  air,  each  having  its  own  equation  of  state.  Each 
region  is  further  divided  into  zones  and  mass  points,  containing  one  half  the 
mass  of  two  adjacent  zones  are  assumed  to  be  at  the  interface  of  these  zones. 


2.1.1  Artificial  Viscosity  Technique 


The  artificial  viscosity  is  a  convenience  first  introduced  by  von  Neumann 
and  Richtmyer  for  the  numerical  treatment  of  shock  waves.  Its  effect  is  similar 
to  the  effect  of  a  real  viscosity.  It  spreads  a  shock  over  a  specified  number 
of  zones  and  thus  permits  one  to  avoid  the  treatment  of  discontinuities  running 
through  discrete  mass  points,  which  presents  serious  difficulties  in  a  finite- 
difference  scheme.  The  spread  of  the  shock  can  be  chosen  and  held  small  by  the  use 
of '  proper  constants.  With  the  constants  used  in  the  present  calculations,  the 
shock  is  actually  spread  over  about  3  zones.  The  artificial  viscosity  is  res¬ 
tricted  to  a  region  which  is  being  compressed  and  is  zero  elsewhere.  The 
original  form  proposed  by  von  Neumann  and  Richtnyer  was: 
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or  its  equivalent  form 


PQ  =  initial  density 
V  =  specific  volume 
Ax  =  Zone  width 
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which  is  valid  only  for  the  planar  case,  where,  u  is  the  flow  velocity. 

6 


Brode  proposed  a  more  general  farm: 
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where,  6=1  for  a  planer  case,  6=2  for  a  cylindrical  case  and  6=3  for  a 
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spherical  case;  R  is  the  distance  from  the  origin:/,  and  A  m  is  the  mass  (per 
steradian,  mass/4ir  for  6  =  3;  mass  per  radian  unit  length,  mass/2ir  i  for 
6  =  2;  mass  per  unit  area,  mass//2  for  6  =  l).  Br ode's  relation  depends  on 
the  geometry.  Wilkins  7,  on  the  other  hand,  gives  the  farm  corresponding  to 
6  =  1  as  valid  for  the  three  geometries .  This  simpler  form  vas  used  since  a 
different  form  of  artificial  viscosity  in  the  chamber  and  in  the  Barrel  might 
lead  to  difficulties,  as  one  zone  can  be  half  in  the  chamber  and  half  in  the 
barrel. 


However,  the  use  of  the  artificial  viscosity  techniques  may  be  strongly 
criticized  as  it  is  applied  in  the  region  near  the  origin.  In  a  spherical 
geometry,  this  technique  is  valid  only  as  long  as  the  shock  thickness  is  small 
compared  with  the  radius  of  the  shock  and  this  is  not  the  case  near  the  origin. 
From  numerous  check  calculations  (Figs.  3  and  4)  it  is  shown  that  in  the  zone 
located  at  the  origin,  the  pressure  due  to  artificial  viscosity  (Q)  takes  on 
values  much  larger  (up  to  about  10  times  in  some  cases)  than  the  flow  pressure 
calculated  from  the  equation  of  state. 


Normally,  in  a  planer  moving  shock  wave  Q  will  be  of  the  order  of  P.  In 
the  center  of  the  shock  near  the  origin  an  additional  compression  due  to  the 
geometry  change  and  the  process  of  reflection  takes  place,  and  consequently  Q 
will  take  on  much  higher  values.  It  is  in  this  region  where  the  total  pres¬ 
sure  (P  +  Q)  is  largely  due  to  the  contribution  of  the  artificial  viscosity, 
that  doubts  arise  whether  the  computation  by  using  (P  +  Q)  will  give  a  reason¬ 
able  approximation  of  the  correct  pressure  history  that  will  act  as  the  driving 
pressure  for  the  shock  or  on  the  projectible  inside  the  barrel.  A  more  detail¬ 
ed  study  of  the  ^  method  near  the  origin  would  be  useful  and  should  be  done. 

This  does  not  seem  to  affect  the  x-t  diagram  of  the  reflection  of  the  im¬ 
ploding  shock,  but  in  all  cases  leads  to  erroneous  values  of  the  pressure  at 
the  origin  and  consequently  it  will  have  some  influence  on  the  Mach  nuaber  of 
the  shock,  or  on  the  velocity  of  the  projectile  in  the  barrel.  For  example, 
from  Fig. 3  at  x  =  5  cm,  and  t  =  45.29  p  sec,  P  =  200  bars  whereas  P  +  Q  =  686 
bars,  that  is,  the  artificial  viscosity  pressure  is  over  twice  the  actual 
pressure.  Consequently,  the  equations  of  motion  will  be  affected.  The  de¬ 
tailed  plot  shown  in  Fig. 3,  indicates  that  in  the  range  50.5  <  t  <  58.0  p  sec, 

Q  >  P.  Consequently,  errors  proportional  to  the  magnitude  of"’ %  compared  with 
the  true  driving  pressure,  P,  can  be  expected. 

2.1.2  Mathematical  Formulation  of  the  Problem 


This  problem  can  be  described  in  a  Lagrangian  form  by  a  set  of  non-linear 
partial  differential  equations  expressing  the  equations  of  mass,  momentum  and 
energy,  as  follows  (Ref .6): 
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Energy: 
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E  =  E(P,V)  for  gas  mixture  2H2  +  0^ 


Equations  of  state :<  P  =  P(E,V)  for  explosive  PETN 

I  E  =  E(P,V)  for  air 


The  artificial- viscosity  pressure,  is  restricted  to  compression  only  and  is 
expressed  as: 
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where  C  is  a  constant  which  determines  the  spreading  of  the  shock  and  is  take*', 
as  C  =  1  for  gases,  hydrogen- oxygen  mixture  and  air,  and  C  =  *s/l.5  for  explosive 
FETN,  which  corresponds  approximately  to  a  spreading  of  3  zones . 


This  set  of  non-linear,  partial  differential  equations  is  then  transformed 
into  a  set  of  finite  differences  equations.  The  two  independent  variables, 
time  and  distance  are  represented  by  N  and  J,  corresponding  to  the  number  of 
the  cycle  being  calculated  and  the  label  assigned  to  a  mass  point.  The  equations 
become: 

Conservation  of  mass  is  assured  since  we  chose  a  fixed  number  of  mass 
points . 
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Ar‘  represents  the  area  of  a  particular  interface,  where  is  the  acce¬ 

leration  of  that  interface.  The  constant  C1  in  the  artificial  viscosity  is 
labelled  CQSQX  in  the  computer  program  (Appendix' 3)  and  is  related  to  C,  the 
usual  constant,  as  C1  =  Q./k  and  mj  is  the  mass  of  zone  J.  We  note  further 
that  the  specific  volume  V  in  the  computer  program  is  normalized  with  respect 
to  the  initial  density. 

It  should  be  noted  that  this  condition  based  on  U  to  determine  if  a  zone 
is  being  compressed  is  only  valid  in  a  onefedimensionai  calculation  and  should 
be  replaced  by  a  condition  based  on  the  specific  volume  which  is  valid  in  all 
geometries  (see  equation  3) •  The  solution  is  then  obtained  by  a  stepwise 
progression  in  time  from  the  acceleration  of  an  interface  based  on  the  old  value 
of  the  pressure.  The  new  velocity  and  position  of  the  interface  can  be  cal¬ 
culated.  A  new  specific  volume  can  then  be  calculated  and  finally  a  simul¬ 
taneous  solution  of  the  energy  equation  and  the  equation  of  state  yields  the 
new  pressure  and  energy.  The  calculation  can  now  be  repeated  for  the  next  time 
step,  etc.  Appendix  B  gives  a  complete  listing  and. a  flow  chart  of  this  code. 

2.1.3  Stability  Requirements 

The  finite  differences  schemes  are  subject  to  mathematical  restrictions 
which  limit  the  size  of  time  increments  that  can  be  taken  without  the  appear¬ 
ance  of  instabilities.  In  the  actual  problem,  we  have  two  stability  conditions: 

a)  C our ant  Condition:  This  condition,  which  is  very  general,  states  only 

that  the  time  increment  At  must  be  small  enough  that  sound  signals  from 
one  mass  point  will  not  have  time  to  reach  the  next  mass  point  during 
this  time  step,  . 

At  <  r  (14) 

where,  c  is  the  speed  of  sound  in  the  zone  and  Ax  its  length.  Appendix 
A  gives  the  detailed  relations  used  for  the  calculation  of  the  speed  of 
sound  in  the  various  regions , 

b)  Artificial  Viscosity  Condition:  In  tne  region  of  compression,  where  the 
artificial-viscosity  pressure  has  a  value,  the  form  of  the  differential 
equations  changes  from  a  wave  equation  to  a  diffusion  type  equation.  It 
is  shown  in  Ref ,6  that  two  successive  time  steps  must  obey  the  following 
relation: 
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2.2  Problem  of  Zoning: 

The  problem  of  zoning  was  previously  found  ^  as  one  of  the  most  important 
difficulties  of  this  code  and  has  not  yet  been  solved  completely  satisfactorily, 
as  a  compromise  is  necessary  between  the  length  of  computing  time,  and  the 
precision  required.  This  problem  involves  in  fact  a  lot  of  other  side  problems: 
the  transition  between  the  hemispherical  geometry  of  the  chamber  and  the  planar 
geometry  of  the  barrel;  the  expansion  of  the  zones  as  they  enter  the  barrel 
resulting  in  extremely  long  zones;  the  computing  time  which  already  approaches 
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10  minutes  on  the  IBM  7094  computer  in  most  calculations  for  the  shock  tube 
case-  Any  refinement  in  zoning  will  increase  this  time  considerably. 

2.2.1  Initial  Zoning 

We  shall  first  review  the  requirements  for  proper  zoning.  Proper  zoning 
must  lead  to  a  solution  of  the  finite  differences  equations  which  is  stable 
and  without  oscillations  in  the  values  of  the  different  physical  variables. 

This  solution  must  be  independent  of  the  number  of  zones  used  in  the  compu¬ 
tation  or  at  least  lead  to  asymptotic  values  as  the  number  of  zones  is  increased 
This  asymptotic  solution  should  be  obtained  with  a  reasonable  number  of  zones 
in  order  to  limit  the  computing  time.  Finally,  this  solution  should  provide 
a  proper  description  of  the  phenomena,  particularly  of  the  flow  in  the  barrel. 
This  effect  on  the  propagating  of  a  shock  has  been  studied  recently  with  this 
code  by  changing  the  initial  zoning  which  will  help  us  to  understand  how  to 
choose  a  proper  zoning  scheme. 

As  a  shock  meets  a  noticeable  increase  in  zone  mass,  some  peculiar  effects 
occur  (see  for  example  the  oscillations  in  the  trajectory  at  t  =  13  nsec  and 
t  =  24  nsec  in  Fig. 5).  The  heavier  mass  point  is  slow  to  accelerate  and  this 
causes  an  increase  in  the  pressure  of  the  zone  located  just  behind  it.  This 
high  pressure  now  accelerates  it  too  fast,  thus  causing  a  decrease  in  the 
pressure  of  this  zone  and  a  decrease  in  the  velocity  of  the  mass  point. 

Finally  the  trajectory  of  the  mass  point  is  oscillating  around  an  average  linr 
and  the  same  happens  to  the  shock.  For  an  important  increase  in  zone  mass, 
the  shock  even  reflects  clearly  on  the  heavier  zone;  a  weak  shock  is  reflected 
and  a  strong  one  is  transmitted.  If  this  zone  mass  were  still  increased 
further,  the  reflected  shock  would  become  the  strongest,  tire  transmitted  shock 
would  become  smaller  and  even,  as  a  limiting  case,  disappear  for  an  infinite 
mass  which  is  equivalent  to  a  solid  wall.  If  we  hrd  a  large  decrease  in  zone 
mass,  a  strong  shock  would  be  transmitted  and  a  rare-f  action  wave  reflected. 

As  the  imploding  shock  wave  reflects  from  the  origin,  these  phenomena  occur 
and  they  cause  the  oscillations  in  the  calculated  pressure.  In  analyzing 
the  different  possible  types  of  zoning,  we  shall  keep  in  mind  these  requirements 
and  compare  how  well  these  zoning  schemes  fulfil  them. 

a)  Constant  Width  Zoning;  In  a  planar,  one-dimensional  code,  constant  width 
zones,  which  are  in  that  case  equivalent  to  constant  mass  zones  are  generally 
used  with  success.  The  first  version  of  this  cod',  which  was  adapted  by 
Piacesi  ®  from  a  planer  code,  also  divided  the  various  regions  into  equal 
width  zones.  This  gave  large  oscillations  in  pressure  at  the  time  of  implosion, 
as  well  as  unacceptable  differences  in  performance  as  various  numbers  of  zones 
were  used  *.  A  constant  width  zoning  leads  to  large  differences  between  the 
mass  of  adjacent  zones  and  this  is  the  reason  for  the  oscillations  which  app¬ 
ear  in  the  pressure  profile  around  the  implosion.  The  constant  width  zoning 
where  the  mass  ratio  of  the  two  zones  nearest  to  the  origin  is  equal  to  8 

was  then  abandoned  and  replaced  by  an  equal  mass  zoning. 

b)  Constant  Mass  Zoning:  This  type  of  zoning  is  without  any  doubt  the  optimum 
one,  relative  to  the  disturbances.  It  gives  no  oscillations  in  pressure  pro¬ 
files,  no  spurious  shocks  reflected  from  a  change  in  zone  mass.  But,  due  to 
the  size  of  the  zone  nearest  to  the  origin  (now  much  larger  if  the  same  total 
number  of  zones  is  used:  400  times  heavier  if  20  zones  are  usrd),  it  results 
in  a  lack  of  detail  as  to  what  is  happening  in  the  barrel  and  consequently 
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results  in  a  dispersion  of  the  final  velocity  of  the  projectile  for  various 
numbers  of  zones.  There  were  then  two  ways  to  overcome  that  difficulty,  either 
to  increase  the  total  number  of  zones  at  the  expense  of  the  computing  time,  or 
try  another  type  of  zoning .  To  increase  the  number  of  zones  was  judged  im¬ 
practical,  since  the  computing  time  is  almost  proportional  to  the  square  of 
the  number  of  zones. 

c)  Sevray's  Zoning:  Sevray  devised  a  new  zoning  scheme  in  which  he  tried 
to  keep  the  main  features  of  a  constant  mass  zoning.  He  divided  the  gas 
regions  in  three  subregions,  each  being  divided  into  constant  mass  zones.  The 
subregion  nearest  to  the  origin  has  smaller  mass  zones.  This  zoning  provides 
sufficient  detail  on  gas  outflow  into  the  barrel  and  consequently  reduces  to 

a  reasonable  value  the  dispersion  in  the  final  velocity  results  as  different 
numbers  of  zones  are  used  (Fig  .6) .  This  zoning  has  been  used  in  the  calculation 
of  the  performance  of  the  UTIAS  hypervelocity  launcher  by  Sevray  ,  and  in  most 
of  the  calculations  of  this  report.  It  is  seen  from  Fig. 6  that  although  the 
details  of  the  projectile  velocity  in  the  barrel  differ  as  the  nuntoer  of  zones 
are  increased, the  final  muzzle  velocity  at  x  =  l60  cm,  is  close  to  14  Km/sec. 

d)  Flagg's  Zoning:  Later,  Flagg  S^consiaering  that  the  differences  in  mass 
from  one  subregion  to  the  next  one  would  introduce  spurious  shocks  and  ex¬ 
pansions,  developed  a  new  partition  of  zone  masses,  which  is  based  on  the 
relation 

^ . +  -  (16) 

i^  +  (n-i)^ 

i  =  1 


where  b  is  the  width  of  the  gas  region,  Ax.  the  width  of  the  jth  zone  and  n 
the  nuacer  of  zones  in  this  region.  ^ 

This  type  of  division  first  appeared  to  have  many  advantages  and  was 
then  used  in  the  calculations  of  the  30  inch  diameter  launched  chamber.  With 
the  same  total  number  of  zones,  the  mass  of  the  zone  nearest  to  the  origin  is 
approximately  10  times  smaller  than  the  mass  of  the  same  zone  in  Sevray's 
calculation,  but  still  10  times  bigger  than  in  a  constant  width  division.  It 
is  thus  providing  a  very  reasonable  detail  of  the  flow  in  the  barrel  without 
slowing  down  the  calculations  too  much.  However,  recent  investigations  have 
shown  oscillations  in  the  pressure  profiles  as  well  as  unreasonably  large 
differences  from  the  results  obtained  by  Sevray,  particularly  in  the  time 
scale  of  the  events.  These  oscillations,  important  mainly  at  the  time  of  the 
implosion  occur  once  again  because  adjacent  zones  have  very  large  difference 
in  mass.  With  the  zoning  of  Sevray,  if  we  have  20  zones  in  the  gas,  the  5 
zones  near  the  origin  are  of  the  oame  mass;  with  the  same  number  of  zones  in 
Flagg's  case,  these  zones  vary  relatively  to  the  corresponding  mass  in  Sevray  s 
case  as  0.10,  0.63,  1-08,  1.55  and  1-85  (a  factor  of  6  between  the  two  first 
zones,  a  factor  of  18  ir.  the  5  first  zones).  With  this  zoning,  as  the  nuntoer 
of  zones  is  increased,  the  time  of  implosion  comes  near  to  the  value  found  in 
Sevray's  calculations.  (That  is  the  time  to  implosion  is  increased  owing  to 
the  cumulative  propagation  time  between  zones).  However,  at  the  same  time  the 
instabilities  m  the  pressure  profiles  are  very  much  increased,  as  shown  on 
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Fig. 7  where  the  same  case  (30  inch  diameter  chamber,  11.65  g  projectile, 

25  kg  of  PETN  and  200  psi  2H2  +  O2  )  is  represented  with  different  zoning 
schemes.  This  increase  in  instabilities  can  be  easily  explained.  The  mass 
ratio  of  the  two  zones  closest  to  the  origin  can  be  expressed  as 


M  -1 
n 

M 

n 


v  -1 
n 


Ax  +  Ax  -1 

_ n _ n 

Sx 

n 


(n-l)' 

- j. 


+  1 


3 


(17) 

An  examination  of  this  rutio  shows  that  it  is  an  increasing  functiomof  r.  for 
n  >  2  varying  between  ) .9?  os  n  =  2,  and  8  as  n  goes  to  infinity. 

Despite  the  advantages  of  this  zoning,  to  obtain  a  detailed  description 
of  the  outflow  in  the  barrel,  this  type  of  zoning  has  to  be  abandoned  as  we 
do  not  know  the  influence  of  the  oscillations  in  the  pressure  at  the  origin 
on  the  final  velocity  and  cannot  tolerate  such  variations  on  the  time  of  the 
implosion. 

Another  zoning  scheme,  similar  in  principle  to  the  one  used  by  Sevray, 
has  also  been  developed  in  order  to  have  smaller  zones  around  the  origin. 

The  gas  region  is  now  divided  into  7  subregions,  each  being  divided  into  the 
same  nuober  of  zones.  Between  each  subregion  there  is  a  mass  ratio  of  the 
order  of  2.  This  permits  one  to  have  zones  around  the  origin  of  the  same 
size  as  the  central  zone  in  Flagg's  zoning,  but  the  fact  that  the  change  in 
mass  is  smaller  and  also  far  from  the  origin  prevents  the  oscillations.  This 
zoning  scheme  will  be  used  further  to  show  the  influence  of  the  mass  located 
at  the  origin  on  the  velocity  profile  of  a  shock  wave  propagating  along  the 
barrel. 

2.2.2  Transition  Between  Spherical  and  Planar  Geometry 

The  treatment  of  the  origin  region,  i.e.,  of  the  transition  between 
the  spherical  geometry  existing  in  the  chamber  and  the  planar  geometry  in 
the  barrel  becomes  a  more  and  more  important  factor  as  the  zoning  scheme  used 
give  smaller  zones  arcurd  the  origin.  When  the  zoning  of  Sevray  is  used,  we 
hacve  only  one  or  two  zones  located  at  the  same  time  in  the  transition  region; 
with  the  finer  zoning,  previously  described,  we  may  have  up  to  about  10  zones 
in  this  region.  More  care  must  then  be  taken  how  to  realize  a  smooth  transi¬ 
tion  from  the  spherical  to  the  planar  flow  and  the  approximate  treatment  of 
Sevray  is  no  longer  acceptable  in  such  a  case. 

It  is  important  to  note  that  this  transition  is  clearly  a  two-dimensional 
phenomenon  and  a  more  complex  type  of  code  should  really  be  used  and  there  is 
no  way  to  treat  it  rigorously  with  a  one-dimensional  code.  Some  approximations 
have  to  be  made  in  order  to  realize  this  transition  in  a  manner  which  does  not 
seriously  affect  the  flow. 

In  the  first  version  of  the  code,  the  flow  is  considered  as  hemispheri¬ 
cal  upstream  from  the  origin  and  planar  downstream.  This  led  to  many  diffi¬ 
culties  and  disturbances  as  shown  previously  in  Ref. 2  and  was  then  abandoned. 
Sevray  made  this  transition  in  two  discontinuous  steps  for  the  area  (Fig. 8) 
but  kept  the  old  routine  for  the  calculation  of  the  specific  volume.  This  is 
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without  consequence  as  ±ong  as  only  one  or  two  zones  are  in  the  transition  region 
tut  if  core  are  inside,  it  must  be  avoided  and  a  smooth  transition  must  be  ad¬ 
opted. 

Later  Flagg  ^  developed  a  continuous  transition ,  taking  the  area  as  the 
one  of  a  spherical  segment  supported  by  the  ba^el  circumference  and  the  volumes 
as  defined  by  these  surfaces.  This  transition  treatment  is  much  superior  to 
Eevray's  treatment,  for  it  is  continuous,  but  it  must  be  noted  that  many  other 
continuous  transitions  could  have  been  developed  and  used  as  well.  All  these 
solutions  have  no  real  physical  meaning  and  are  somewhat  arbitrary,  but  the 
treatment  of  a  two-dimensional  phenomenon  with  a  one-dimensional  code  can  be 
nothing  but  approximate  ar.d  only  a  two-dimensional  code  can  treat  it  rigorously. 
This  problem  is  presently  being  considered  at  UTIAS. 

2.2. 3  Redivision  of  Long  Zones  in  the  Barrel 

The  barrel  is  of  a  very  small  diameter,  so  that  for  any  reasonable  number 
of  zones,  one  zone  as  it  expands  can  occupy  a  very  long  distance  in  the  barrel, 
thus  providing  insufficient  detail  on  what  is  happening  in  the  barrel.  To 
overcome  this  problem,  Sevray  introduced  a  scheme  to  split  the  zone  next  to  the 
projectile  in  two  zones  of  equal  mass  and  pressure,  whenever  its  length  was 
greater  than  3  cm.  This  scheme  was  approximate  and  worked  very  well  as  long  as 
a  projectile  case  was  calculated.  Later,  as  the  shock-tube  cases  were  run,  it 
led  to  very  big  disturbances  and  another  scheme  had  to  be  found. 

To  explain  why  this  happened,  let  us  consider  the  momentum  equation  for 
the  projectile* 


In  the  denominator  of  this  expression,  m  ^  ,  is  the  most  important  term 
and  the  division  by  2  of  n .  i_,  as  the  zone  is  dfyiaed,  has  no  serious  influence 
on  the  acceleration  of  the* projectile.  However,  if  we  now  consider  a  shock 
tube  case  and  look  again  at  the  denominator  of  this  expression  applied  at  the 
interface,  we  note  that,  nipr0j  being  equal  to  zero  and  m.+i  being  very  small, 
owing  to  the  low  pressure  in  the  channel,  m^+1  becomes  the  most  important  term, 
so  that  dividing  i.  by  2  is  almost  equivalent  to  multiplying  the  acceleration 
of  this  interface  oy^2,  thus  leading  to  the  disturbances  we  found.  Even  l 
these  disturbarces  damp  out,  what  their  actual  influence  on  the  result  maybe 
is  difficult  to  evaluate  and  it  is  preferable  to  avoid  them. 


Another  scheme  based  on  a  more  realistic  hypothesis  was  developed.  It 
permits  a  splJt  of  any  zone  in  the  barrel  into  two  zones  of  equal  mass  (and  not 
only  the  zone  placed  just  behind  the  interface)  whenever  its  length  increases 
toe  much.  Essentially,  this  scheme  keeps  constant  the  values  of  the  acceleration 
cn  both  boundaries  of  the  zone  which  will  be  divided,  and  that  permits  one  to 
assign  to  each  new  zone  a  pressure  calculated  from  the  momentum  equation.  The 
specific  volumes,  rather  than  being  taken  equal  in  both  zones,  are  now  calculated 
from  the  pressure  assuming  a  y  law  relation  between  pressure  and  density  and  an 
effective  y  equal  to  1.1b  as  mentioned  by  Flagg  9 ,  The  energy  is  then  calculated 
from  pressure  ard  density  using  the  equation  of  state,  and  the  conservation  cf 
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energy  gives  the  kinetic  energy  of  the  new  mass  point  and  thus  its  velocity. 

This  procedure  permitted  one  to  obtain  a  much  better  description  of  the  flow  in 
the  barrel,  but  in  consideration  of  the  other  limitations  of  the  code  and  of 
the  increase  in  computing  time,  it  was  used^only  in  a  small  number  of  cases 
to  check  its  influence.  This  influence  was  found  to  be  not  important  and  justi¬ 
fies  the  fact  that  it  was  not  used  regularly  in  the  calculations. 

2.3  Equations  of  State 

Another  major  difficulty  of  such  calculations  is  to  find  appropriate 
equations  of  state  for  the  different  materials  involved,  particularly  when  the 
conditions  are  extreme  as  they  are  in  our  case  where  the  pressure  and  temperature 
at  the  origin  are  ideally  infinite  at  the  time  of  implosion.  The  three  equations 
of  state  that  were,  used  are  real  gas  equations  of  state,  but  have  been  used  well 
beyond  their  range  of  applicability. 

2.3.1  Equation  of  State  of  the  Gas  Mixture  2  v2 

It  is  an  extrapolation  of  the  real  gas  calculation  of  Moffatt  ^  done  by 
Brode  and  programmed  by  Piacesi.  It  is  expressed  as 

E  =  6.57  x  PV  +  _  974.0  (IV)  _  a  |0.101xlo-3Ln  _P -  .  o.2325xlO~ 3 

1140.0  +  (PV)  *-  1.013xl03 

(19) 

where: 

a  =  0  for  PV  <  1.0465 

a  =  8600  PV-9000  for  1.0465<PV  <  3.488 
a  =  21xl03  for  3.488  <  PV 


the  units  are: 

P  in  10 10  erg/cc 
E  in  1010  erg/g 
V  in  cc/g. 

This  equation  is  representative  of  the  21^+  02  mixture  in  the  range: 

l600°  K  <  T  <  6000°  K 
0.01  bar  <  P  <  1000  bars 

and  has  been  used  well  beyond  these  limits  in  our  calculations,  since  unfort¬ 
unately  no  calculations  are  available  for  a  stoichiometric  mixture  of  hydrogen- 
oxygen  for  temperatures  to  10°  °K  and  pressure  of  the  order  of  Kr  bars.  How¬ 
ever,  it  has  been  shown  that  at  high  pressure  this  equation  of  state  is  closely 
equivalent  to  a  flow  with  7  =  I.Ih,  which  is  in  good  agreement  with  other 
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=  10  bars 
=  10*  bar-cc/g. 
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available  calculations 


2.3.2  Equation  of  State  for  the  Explosive 


This  equation  is  a  numerical  fit  to  the  experimental  data  of  a  50^ 
mixture  of  PEIN  and  TNT  (pentolite)  provided  and  programmed  by  Piacesi.  It  is 
written  as 


P  (E,  V)  =  A 


(p)  *  B(p)  E 


(20) 


where 


A  =  0.002l64p/4  -  2.0755e  *  6*°/P 


B  =  0.35P 
P. 


o 

p  *  “ 


po  is  the  initial  density  and  was  taken  as  PQ  =  O.58  g/cc.  The  units  are: 


P  in  megabars 


E  in  megabar  - 


-c/& 


It  must  be  mentioned  that  both  the  equation  of  state  for  the  gas  and 
the  equation  of  state  for  the  explosive  do  not  give  directly  the  temperature  and 
in  this  code  we  deduce  it  from  the  ideal  relation  FV  =  RT  where,  R  =Cl/M  and 
M,  the  molecular  weight.  The  molecular  weight  of  the  mixture  varies  with  the 
degrees  of  dissociation  and  ionization  and  is  thus  a  function  of  the  temperature . 
For  the  2Hg+  02  mixture,  M  was  taken  as  12  which  corresponds  to  the  unreacted 
mixture.  For  the  explosive,  M  was  taken  as  1,  since  we  did  not  know  the 
composition  of  the  products  of  detonation  of  pentolite  and  corrections  have  to 
be  made  to  obtain  a  reasonable  value  of  the  temperature.  The  temperatures  are 
thus  inaccurate  and  give  only  some  indications  of  their  variations  with  time 
or  position.  It  is  hoped  that  better  calculated  temperatures  will  be  obtained 
in  future  calculations,  at  IZTIAS. 


2.3.3  Equation  of  State  for  Air 


First  a  perfect  gas  equation  FV  =  RT  was  used,  but  later  it  was  found 
a  more  accurate  equation  would  be  preferable.  An  equation  of  state  was  then 


programmed  from  an  extrapolation  of  the  calculations  of  Gilmore  ^  done  by 
Brode  ^-3. 


A  convenient  form  to  write  this  equation  of  state  for  air  is 

e  =  W) 


7-1 

e  =  5  (m-1) 


(21) 


M  =  m(C  »Tt) 
|  =  Ln  (t,) 
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with  the  non-dimensionaiized  variables 

p 

tt  =  —  with  P  =  1  atm. 

Po 

n  =  £  with  Po  =  1.293  10"3  g/cc. 
Ho 

E 

e  =  —  with  E  =  i960.  0  bar-cc/g. 
o 

6  =  |  With  To  =  273.92  °K. 

o 


we  have 

S(C) 


485  +  3860 

1000  +  C  7500  +  16.5C 


(22) 


h  =  hQ  +  0.09  (h0-h2)  i  (23) 

taking  y  =  l/5,  and  are  given  by  the  following  relations 

w  .  j  +  25.894868  y+j  +  §60  2536  y(l:y)  +  jjiopoyjky).  (pt) 

Mo  4.778973  y+1  3000  yz+l  9  104  y*  +1  12  10b  y2  +1  1 ’  ' 


6002  y  +  4 
1000  y+1 


(25) 


this  equation  of  state  is  valid  for  temperatures  up  to  24,000  °K  and  was 
used  almost  only  in  its  range  of  validity. 

2.4  Detonation  Scheme 

In  programming  the  equation  s  of  state  for  explosive  and  for  gas  mixture 
2H2  +  02  the  detonation  scheme  described  by  Wilkins  ^  was  used.  In  each  zone, 
the  burn  fraction  F  defined  as 


is  calculated  and  the  pressure  in  the  zone  is  defined  as 

P  =  P  (E,  V)  x  F  (27) 

where,  V  is  the  specific  volume  and  Vcj  the  Chapman-Jouguet  volume.  This 
burn  fraction  F  is  used  to  spread  the  detonation  front  over  several  zones  in 
the  same  manner  as  the  artificial  viscosity  spreads  a  shock  front  over  several 
zones.  The  burn  calculation  is  started  by  setting  F  =  1  in  the  zone  that 
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corresponds  to  the  point  of  initiation  of  the  detonation.  When  F  =  i  for 
a  zone,  then  ail  the  chemical  energy  contained  in  that  zone  has  been  deposited 
and  F  remains  equal  to  1  thereafter. 

It  should  be  noted  that  the  calculations  will  proceed  over  several 
zones  (Wilkins  '  mentions  2  to  3  times  the  number  of  zones  on  which  the 
artificial  viscosity  is  spread  over,  i.e.,  15  zones  approximately  in  our  prob¬ 
lem,  but  calculations  seem  to  show  a  much  smaller  number)  before  the  deton¬ 
ation  front  is  correctly  established.  This  implies  the  use  of  a  maximum 
number  of  zones  in  order  to  obtain  reasonable  accuracy  in  the  detonation  velo¬ 
city  and  consequently  in  the  time  scale  of  the  phenomenon.  The  Chapman- 
Jouguet  volumes  which  were  used  in  this  work  are 

V  .  (PETN)  =  0.7872  (28) 

C  J 

VcJ  (2H2  +  02)  =  0.54 

p_-3 

In  previous  works  J,  unfortunately  an  erroneous  value  of  the  Chapman- 
Jouguet  volume  for  PETN  has  been  used.  Some  runs  have  been  done  to  check  the 
influence  of  this  parameter  and  it  was  found  this  error  did  not  affect  the 
results  seriously  as  shown  in  Figs.  9  and  10.  It  is  seen  that  the  correct 
value  (VCJ  =  O.7872)  gives  a  higher  (4<£)  final  contact  surface  velocity  and 
lower  implosion  pressure  at  the  origin. 

3*  RESULTS  OF  THE  COMPUTATIONS: 

Only  a  few  results  indicative  of  the  general  trend  will  be  given  in  this 
report  since  it  has  been  recently  found  follcwing  a  comparison  with  experi¬ 
mental  results,  that  they  were  of  limited  valu».  In  Section  5>  the  reasons 
for  this  fact  will  be  explained  in  detail  and  some  indications  of  the  remedy 
that  could  be  applied  to  improve  these  calculations  and  make  them  more  re¬ 
liable  will  be  given. 

3-1  The  Shock  Tube  as  a  Limiting  Case  of  a  Hypervelocity  Launcher  as  the 
Mass  of  the  Projectile  Goes  to  Zero 

The  case  presented  consists  of  an  explosive  loading  of  100  g  PETN,  an 
initial  pressure  in  the  chamber  of  hOO  PSIA  2Hg  +  O2  and  an  air  pressure  of  1 
torr  in  the  barrel.  Four  runs  were  made  with  the  mass  of  the  projectile 
varying:  1.72  g,  1.0  g,  0.5g  and  0.2  g.  The  results  are  shown  on  Figures  11 
and  ]2.  Figure  II  indicates  the  variation  of  the  velocity  of  the  projectile 
or  interface  at  2  meters  as  the  mass  of  the  projectile  decreases  and  goes  to 
zero.  On  Figure  12  is  giver,  the  variation  of  the  velocity  profile  at  a  dis¬ 
tance  of  two  meters  from  the  origin  with  the  mass  of  the  projectile  as  a 
parameter. 

From  Figure  11  it  can  be  seen  that  as  the  mass  of  the  projectile  de¬ 
creases,  the  importance  of  the  first  peak  pressure  on  the  base  of  the  pro¬ 
jectile  becomes  bigger  and  finally  as  the  mass  decreases  still  more,  the 
influence  of  the  second  peak  disappears  completely.  The  trajectory  of  the 
projectile  is  now  determined  by  the  first  implosion  and  shock  waves  in  the 
barrel  generated  by  the  following  implosions  will  no  longer  overtake  it. 

Figure  12  shows  that  as  expected  the  projectile  velocity  increases  with  de¬ 
creasing  mass.  It  was  not  possible  to  go  to  masses  m  <  0.2  g,  as  the  program 
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-  c u.:-!C_e  give  a  -.ou’ticn  :o 
-ou-.cl  yield  the  chock  tube  case, 
presently  urnier  review. 

3.2  dlrvK.  Tube  Calculation 
;,.vM  Gdc  Case: 

Results  were  obtained  for  gas  cases  with  an  initial  pressure  in  the 
cheooer  of  2CC,  400,  Co 0,  80C  PSL4  2ru  +  (V  .  Figure  13  gives  the  different 
■velocity  profiles  obtained.  It  is  seen  they  are  very  similar  and  the  velo¬ 
cities  of  the  interfaces  differ  only  slightly.  This  is  not  physically  real¬ 
istic  for  we  can  expect  a  considerable  increase  with  initial  driving  pressure 
especially  near  the  origin,  before  attenuation  sets  in.  The  fact  that  no 
c.-.^nge  is  observed  is  indicative  of  some  computational  problem.  These  un¬ 
expected  computational  results  coupled  with  the  experimental  results  show 
(see  Section  6)  the  inability  of  this  code  to  describe  the  phenomenon. 

3.2.2  g>' plosive  Case: 


r  era:-. Li  masses.  1  hi s  i ;  unforuun.: ie 
The  entire  num/jrica'.  program  is 


A  series  of  runs  with  tne  8  inch  diameter  chamber  and  a  5/16  inch 
diameter  barrel  has  been  run  for  a  constant  initial  pressure  of  200  PS1A  in 
the  chamber  anu  diverse  explosive  loading:  100  g,  200  g  and  4o 0  g.  Figure 
14  gives  the  interface  velocity  for  these  diverse  cases.  It  is  seen  that  the 
final  velocity  at  10  meters  first  increases  with  increasing  explosive  weight, 
goes  through  a  maximum  and  then  decreases  as  the  weight  of  explosive  is 
further  increased.  This  emphasizes  the  need  for  an  optimization  of  the 
initial  conditions  tc  obtain  the  highest  Mach  number  (see  Hex'- 2).  However, 
in  the  light  of  the  comparison  between  computed  and  experimental  results, 
this  optimization  has  not  been  completed  using  the  present  program. 


Figure  15  is  the  tirar— distance  diagram  of  a  shock  tube  case  with  an 
explosive  loading  of  200  g  PETIf  and  an  initial  pressure  of  200  K>IA  2R3+  O2 
for  2  5/ l6  inch  diameter  channel.  The  successive  implosions  were  followed 
and  it  is  interesting  to  remark  that  the  successive  shocks  they  generate  in 
the  barrel  propagate  with  lower  and  lower  velocity  and  will  never  overtake 
sc  that  the  trajectory  of  the  interface  is  in  fact  determined  only  by  the 
first  implosion.  It  is  seen  that  shock  Mach  numbers,  M  ~  100  may  be  achieved. 


Figure  16  gives  the  pressure  profiles  at  various  instants  corresponding 
the  case  shown  on  Fig. 15-  It  is  seen  that  uniform  flow  exists  between  the 
contact  surface  and  the  shock  wave.  If  this  result  is  verified  experimentally 
it  would  be  most  encouraging  t.c  use  the  present  facility  to  generate  uniform 
flow  regions  at  extreme  temperatures  and  pressures.  Figure  17  shows  the 
luluence  of  counterpresxure  on  the  velocity  of  the  interface.  There  is  only 
a  few  per  cent  difference  between  0.1  torr  2nd  1  torr.  However,  the  profile 
consider abl.y  and  inexplicably  as  the  counter pre s sure  is  increased  to 
V-i  torr.  The  remarks  made  in  the  previous  subsection  also  apply  here.  It 
can  te  seen  chat  during  the  early  phase  of  the  flow  the  counterpressure  has 
tittle  effect  on  the  accelerating  interface  velocity.  Whereas,  in  practice, 
one  gar  escape  speed  as  obtained  from  the  chamber  temperature  sets  an  upper 
•Jiund  near  the  origin  and  viscous  effects  will  decay  this  speed  with  distance. 


•  •  KfSRjaaigAL  aarji/nd: 

.  1  Id  iuipment 

- i  1  experimeic 3 1  runs  were  done  in  an  8  inch  diameter. 


hemispherical 


chamber.  The  barrel  was  mace  of  stainless- steel,  high-pressure  tubing  with  a 
5/16  inch  internal  diameter  and  was  approximately  4  meters  long.  Scribed 
stainless  steel  diaphragms,  0.015  inch  thick,  which  had  been  tested  previously 
by  Watson  were  used  and.  chosen  following  his  recommendations. 

The  velocity  of  the  shock  was  measured  using  three  different  methods: 

1)  Ionization  Gauges:,  they  detect  the  arrival  of  the  ionization  front 
which  accompanies  the  shock.  They  were  the  principal  measuring  system. 

2)  Photomultipliers:  they  detect  the  arrival  of  the  luminous  front  which 
accompanies  the  shock.  They  have  been  used  as  an  independent  way  of 
measuring  the  time  of  arrival  of  the  shock  at  a  fixed  station  in  order 
to  check  the  results  given  by  the  ionization  gauge  technique. 

3)  Microwave:  in  association  with  Elsenaar  this  method  should  have 
given  a  continuous  measurement  of  the  velocity  of  the  shock  wave  and 
by  following  the  profile  it  should  have  been  possible  to  evaluate 
boundary- layer  effects  and  shock-wave  attenuation  in  the  barrel.  But 
unfortunately  the  reflection  of  the  microwaves  from  the  shock  was  very 
weak  owing  to  the  low  power  level  of  the  existing  equipment  and  no 
practical  results  could  be  obtained. 

4.2  Ionization  Gauges 

The  ionization  gauges  consisted  of  a  thin  copper  wire  (l/32  inch 
diameter)  placed  in  a  hole  in  the  wall  of  the  barrel  and  well  insulated  from 
the  wall.  They  extended  about  0.5  mm  inside  the  barrel.  Over  each  ionization 
gauge,  the  barrel  was  reinforced  with  a  tight  collar  to  prevent  the  gauge 
from  leaving  its  hole  as  the  pressure  in  the  channel  increased,  following  the 
arrival  of  a  shock.  Five  ionization  gauges  were  placed  along  the  barrel  at 
fixed  intervals  of  14  inches.  During  an  experiment,  they  were  charged  at 
300  volts  D.C.  To  prevent  accidental  discharge,  they  were  checked  before  a 
run  at  500  volts  D.C.  These  ionization  gauges  were  connected  to  a  circuit 
transforming  their  signals  into  sharp  pulses  which  were  recorded  on  an  osci¬ 
llograph.  We  had  already  at  TJTIAS  one  circuit  (Fig.l8)  designed  by  Flagg  9 
for  his  measurement  of  the  incident  detonation  wave  and  reflected  shock  in 
the  one-dimensional  chamber.  This  circuit  had  been  designed  to  be  able  to 
distinguish  shocks  at  intervals  of  about  5  ms.  This  led  to  difficulties  in 
the  first  measurements.  The  same  ionization  gauge  was  giving  several  signals 
corresponding  to  the  arrival  of  the  different  shocks  caused  by  the  successive 
implosions  and  on  the  record  uhere  was  no  way  to  distinguish  between  them. 

A  new  circuit  was  then  designed  with  a  time  constant  for  charging  of  the  order 
of  150  milliseconds  and  a  time  constant  for  discharging  of  the  order  of  2 
microseconds.  This  circuit  represented  on  Fig.l8  was  built  and  used  with 
success.  The  signals  attenuated  28.2  times  were  recorded  on  a  zig-zag  (  or 
raster)  oscillograph  together  with  the  signals  of  a  10-microseconds  time-mark 
generator.  The  oscillograph  was  triggered  from  the  ignition  and  a  time  of 
76  microseconds  was  assumed  between  the  ignition  and  the  first  implosion. 

4 . 3  Verification  of  the  Ionization  Gauge  Technique 

After  the  first  experiments  which  showed  velocity  profiles  quite 
different  from  the  computed  profiles,  it  was  decided  to  check  the  ionization 
gauge  technique.  This  was  done  by  placing  at  the  same  distance  from  the 
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origin  an  ionization  gauge  and  a  photomultiplier  looking  through  a  window  into 
the  barrel.  The  window  was  a  l/l6  inch  diameter  hole,  drilled  perpendicularly 
to  the  barrel  and  filled  with  a  transparent  resin.  The  time  of  arrival  of  the 
luminous  front  and  of  discharge  of  the  ionization  gauge  were  recorded  separately 
as  shown  on  Fig. 19  for  a  200  PSIA  +  O2  case.  The  agreement  between  the  two 
techniques  was  found  to  be  very  good. 

b .4  Experimental  Results 

Only  a  few  cases  have  been  run  so  far,  mainly  for  two  reasons: 

1)  They  showed  sufficient  evidence  that  the  computations  were  not  adequate. 

2)  We  were  strongly  limited  in  our  choice  of  initial  conditions  by  the 
strength  of  the  hemispherical  chamber,  which  had  been  previously 
damaged  and  could  no  longer  stand  explosive  runs. 

For  these  reasons  only  three  cases  have  been  run  and  duplicated  in 
order  to  obtain  average  values  since  the  reproductibility  was  not  always  very 
good;  initial  pressure  of  these  were  200  K3IA,  hOO  PSIA  and  600  PSIA  2H^  ■  0g 
without  explosive  and  with  a  counterpressure  of  1  torr  or  10  torr.  ^ 

Figure  20  gives  the  oscillograph  recordings  of  a  bOO  PSIA  2Hg  +  0g 
case.  On  Fig.  21  the  velocity  profiles  versus  distance  of  the  three  studied 
cases  are  given.  Two  points,  where  the  velocity  of  the  shock  in  the  400  PSIA 
case  was  measured  with  the  microwave  equipment  -*-6  are  also  given  on  this 
previous  figure.  These  points  are  in  good  agreement  with  the  ionization  gauge 
measurements .  The  three  profiles  have  the  same  general  character:  very  high 
initial  velocity  followed  by  an  extremely  rapi  l  attenuation  of  the  shock.  For 
example,  in  the  600  psi  2Hg  +  0 2  case,  the  she.  k  wave  attenuates  from  Ms  ~  40 
to  Ms  ~  10,  in  6  ft.  This  attenuation  was  firs*  judged  too  strong,  but  it  is 
in  fact  not  so  strange  when  we  consider  such  high  Mach  numbers  flow  in  a  pipe 
of  such  a  small  diameter,  and  we  have  probably  a  fully  developed  pipe  flow. 

Other  runs  have  been  done  last  year  by  Flagg,  with  variable  explosive  loading 
since  at  that  time  tht  chamber  could  stand  200  g.  PETN.  These  runs  show  the 
same  general  features  but  they  have  been  done  with  only  two  ionization  gauges 
and  consequently  give  only  a  rough  idea  of  the  velocity  of  the  shock. 

More  experimental  runs  will  have  to  be  done  in  order  to  evaluate 
completely  the  possibilities  of  this  device  and  the  influence  of  radiative, 
convective,  ablative  and  frictional  losses,  and  to  optimize  its  performance. 

With  the  new  launcher  which  has  been  designed,  it  will  be  possible  to  use  a 
1  inch  diameter  barrel  and  this  should  permit  us  in  the  future  to  carry  out 
more  realistic  and  useful  experiments.  However,  for  the  aim  of  this  study, 
these  few  runs  have  been  sufficient,  permitting  us  to  test  well  the  code  used 
and  to  show  its  inability  to  describe  accurately  the  phenomena  occurring  in 
such  a  shock  tube  as  will  be  shown  subsequently. 

5.  COMPARISON  BETWEEN  THEORETICAL  AM)  EXPERIMENTAL  RESULTS 

Figure  22  shows  cn  the  same  plot  the  velocity  profile  of  the  shock  in  a 
200  PSIA  2H2  +  0g  gas  case  computed  from  this  code  and  as  measured  experimentally. 
The  two  profiles  are  quite  different  and  it  was  suspected  that  the  ionization 
gauges  were  not  working  properly.  However,  after  the  gauges  were  checked  by  an 
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independent  technique,  the  optical  measurement  of  the  arrival  of  the  luminous 
shock  front,  and  after  having  duplicated  all  runs  and  obtained  results  within 
20  per  cent  it  was  no  longer  possible  to  doubt  the  results.  The  calculations 
were  then  checked  with  great  care  and  we  arrived  at  the  conclusion  that  the 
computer  code  was  unable  to  describe  accurately  the  starting  process  of  this 
shock  tube. 

Let  us  consider  the  movement  of  the  interface  between  the  gas  mixture 
2H2  +  O2  in  the  chamber  and  the  air  in  the  barrel.  Just  at  the  time  at  which 
the  diaphragm  breaks.  This  movement  is  governed  by  the  momentum  equation: 

A”  [ 
i  £mj.§  + 

For  a  given  pressure  difference  across  the  interface,  the  acceleration  depends 
on  the  values  of  the  masses  M.  1  and  M.+i.  In  the  channel,  the  pressure  is 
very  low,  so  that  the  mass  M^i2  is  very2 small.:  typically,  for  a  pressure  of 
1  torr  and  a  5  cm  long  zone,  2  Mj+i_  is  equal  to  O.h  10“5g.  On  the  other  side, 
M*  1,  the  last  zone  in  the  gas  mixture  is  in  a  typical  run  using  Sevray’s 
zoning  2  with  20  zones  in  the  gas  of  the  order  of  0.1  g,  so  that  the  accele¬ 
ration  of  the  interface  is  almost’  inversely  proportional  to  the  mass  of  this 
zone.  This  picture  is  rigorously  true  only  for  the  cycle  at  which  the  dia¬ 
phragm  will  burst,  as  the  breaking  pressure  has  been  constantly  chosen  as  1000 
bars  .■  For  subsequent  cycles ,  the  driving  pressure  will  also  be  dependent 
on  the  mass  Mj + 1_ .  However,  the  influence  of  the  value  of  the  mass  located  at 
the  origin  was  2clearly  shown  as  another  type  of  zoning  problem  in  Sec.  2.2.1. 
The  zoning  used  allowed  us  to  have  a  mass  at  the  origin  approximately  10  times 
smaller  than  with  Sevray's  zoning,  with  only  35  zones.  Figure  23  shows  how 
the  velocity  profile  varies.  As  expected  from  previous  consideration,  the 
inital  acceleration  of  the  interface  is  much  higher.  In  addition,  the  final 
velocity  is  much  higher.  The  latter  result  could  not  have  been  predicted 
beforehand.  This  can  be  understood  and  explained:  in  the  first  case,  with  a 
heavy  mass  at  the  origin,  it  takes  a  long  time  for  the  interface  to  attain  its 
final  velocity  and  is  still  accelerating  at  10  meters  from  the  origin,  where 
the  calculation  was  stopped;  when  a  much  lighter  mass  is  located  at  the  origin, 
the  interface  accelerates  much  faster  and  the  final  velocity  is  obtained  in 
a  much  smaller  distance.  It  is  interesting  to  mention  that  if  we  apply  the 
theory  of  a  perfect  shock  tube  ^7 ,  the  initial  acceleration  of  the  interface 
is  infinite.  From  the  previous  considerations  it  is  seen  that  to  approach 
such  a  result  we  should  have  Mj_A.  of  the  same  order  as  Mj+i.  However,  the 
continuity  of  the  masses  across  uhe.  interface  is  difficult2to  obtain  and  would 
require  an  extremely  large  number  of  zones. 

The  results  of  Section  3*2.1  can  now  be  explained.  The  shock  Mach 
numbers  calculated  for  initial  pressures  of  200,  400,  600  and  800  PSIA 
2H2  +  O2  were  found  to  be  almost  the  same,  as  shown  on  Fig.  13,  when  they  should 
have  increased  with  initial  pressure.  In  the  computations  when  the  initial 
pressure  is  doubled,  the  pressure  ir.  the  chamber  is  also  approximately  doubled. 
However,  if  we  keep  the  same  total  number  of  zones  in  the  gas,  the  mass  of 
the  zone  located  at  the  origin  is  also  doubled  end  both  effects  compensate  in 
the  momentum  equation  to  give  the  same  acceleration. 

From  the  previous  discussion,  it  is  seen  that, 
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a)  the  computed  velocity  profile  is  not  realistic 

b)  the  velocity  obtained  at  10  meters  is  not  the  final  velocity, 
since  the  interface  is  still  accelerating 

Consequently,  the  code,  as  used  is  inadequate  for  such  calculations.  The  follow¬ 
ing  ir;... movements  may  be  implemented  in  order  to  make  the  calculations  closer  to 
reality. 

6.  FUTURE  V.’ORK 


Severui  changes  to  the  present  code  will  have  to  be  made,  but  there  are 
three  main  directions  in  which  improvements  will  have  to  be  found  and  these  will 
be  reviewed  in  order  of  importance. 

a)  Initial  Zoning:  As  shown  in  Section  2.2.1  and  5,  large  discontinuities  in 
zone  mass  cannon  be  tolerated  and  continuity  of  zone  mass  must  be  insured  across 
the  i  verfaoe.  From  recent  discussions  with  Dr.  R.  L.  Erode  it  has  been 
"‘Oner nded  that  even  Sevray's  zoning  is  not  acceptable  since  the  ratio  of  zone 
masses  of  two  adjacent  subregions  vnich  is  equal  to  3  is  too  big.  However,  a 

lu  per  cent  increase  in  mass  from  cne  zone  to  the  next  may  be  acceptable.  This 
suggests  the  use  in  the  chantoer  of  zoning  in  which  as  we  go  away  from  the  origin, 
the  zone  mass  is  increased  by  a  fixed  percentage  of  the  order  of  10  per  cent. 

A  reasonable  number  of  zones  would  have  to  be  assumed  in  the  explosive  region 
and  the  number  of  zones  in  the  gas  region  would  be  determined  by  requiring  con¬ 
tinuity  of  zone  mass  through  the  interface  gas  mixture  2H2  +  02/air  and  gas 
mixture  2H2  +  02/explosive.  Such  zoning  would  be  ideal.  There  are  no  longer 
any  important  discontinuities  in  zone  mass  and  thus  spurious  shocks  or  expansions 
ere  avoided.  The  mass  located  at  the  origin  is  now  very  small  and  of  the  same 
size  as  the  masses  in  the  barrel  thus  being  independent  of  the  pressure  in  the 
chamber,  and  the  difficulties  explained  at  the  end  of  Section  5  are  avoided. 
However,  such  zoning  will  lead  to  very  large  numbers  of  zones  and  require 
extremely  large  storage  of  data.  The  time  of  computation  will  be  very  long  and 
will  limit  the  number  of  cases  that  can  be  run,  so  that  a  compromise  will  prob¬ 
ably  have  to  be  found.  An  additional  problem  is  that  for  very  small  zones  the 
transition  between  t\;  spherical  and  planar  geometry  (see  Sec. 2. 2. 2)  becomes 
extremely  ltcpor‘.uut.  As  noted  before,  only  a  two-dimensional  treatment  of  the 
transition  region  will  describe  the  local  flow  accurately. 

b)  Equations  of  State:  When  a  more  accurate  code  will  have  been  developed,  it 
will  become  necessary  to  use  more  precise  equations  of  state  in>  order  to  be  able 
to  compare  expe.imental  and  computed  results.  The  present  equations  have  been 
used  well  beyond  t.heii  limits  of  validity.  This  problem  is  actually  under 
investigation  and  it  is  hoped  that  we  will  soon  have  a  better  equation  of  state 
for  the  mixture  2E-,  +  0o.  The  equation  of  state  for  the  explosive,  PETN,  will 
also  have  to  be  reconsidered  since  the  equation  we  use  is  in  fact  an  equation  for 
pentolite  (50  per  cent  of  PETN  and  of  TWT). 

c)  Artificial  Viscosity  Technique:  The  Q,-method  will  have  to  be  investigated 
much  “ore/  The  form  of  the  artificial  viscosity  should  be  chosen  as  described 
by  Brode  °  and  the  criterion  presently  based  on  the  velocities  will  have  to  be 
replaced  by  a  criterion  based  on  the  specific  volume.  The  artificial  viscosity 
technique  leads  to  an  erroneous  value  of  pressure  in  the  zone  located  at  the 
origin  and  this  effect  on  the  final  performance  will  have  to  be  evaluated. 
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CONCLUSIONS 


7. 


The  experimental  measurements  of  shock  velocity  using  ionization  gauge 
have  shown  that  the  numerical  code  vised  in  this  work  was  riot  very  successful 
in  predicting  the  performance  of  the  UTIAS  implosion-driven  shock  tube.  The 
deficiencies  of  the  code  have  been  analyzed  in  detail  and  recommendations  have 
been  given  to  overcome  some  of  these  difficulties.  The  main  problem  is  to 
choose  proper  initial  zoning  that  would  avoid  computational  effects  and  errors 
arising  from  a  large  change  in  zone  mass  and  which  will  still  not  lead  to 
unreasonable  computing  times. 
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FIG.  5  INFLUENCE  OF  A  CHANGE 
IN  ZONE  MASS 


FIG.  9  INFLUENCE  OF  THE  VALUE  OF  THE 
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FIG.  18  IONIZATION  PROBE  CIRCUITS 
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APPENDIX  A:  SOUND  SPEED  RELATIONSHIPS 


2H2+02  mixture 

The  equation  of  state  is  of  the  form 

E  =  E  (P,V) 

The  speed  of  sound  in  a  gas  is  given  by  the  relation 


The  total  derivative  for  P  is 

dP  =  (sv)E  dV+  (li)v  dE 

and  thus: 

(§0.  -(&HHCSO. 

The  thermodynamic  partial  derivatives  are  related  by 

(!*),(!),(«),■•■ 

and 


substituting  we  obtain 

(If  )  =  -p  ('i-)v-  p  \ 
'•[P  )v+  (^i  )v(s  \  1 

-  K!)J/(!)V 


A1 


APPENDIX  3 

PROGRAM  LISTING  AND  FLOW  CHARTS 


TABLE  1 

Summary  of  Subroutines 

DECK 

SUBROUTINE 

PAGE 

NAME 

NAME 

NUMBER 

MAIN 

1 

READ 

LIRE 

2 

SETUP 

DEPOT 

3-6 

GROS 

GUIS 

7-8 

RETOUR 

LFTOV 

9-10 

OUTPUT 

SORT  I 

11-12 

ETAT  1 

EQST  1 

13-14 

ETAT  2 

EQST  2 

15 

ETAT  3 

EQST  3 

16 

APPENDIX  B 


PROGRAM  LISTING  AND  FLOW  CHARTS 


SIAFTC  MAIM  NOLIST 

COMMON  A^EAl  .PRAKE.CSQVA  »DLAVA .DTLAS  .D'JO.El  NSU.EKSUM .EMPRO.ESUM* 

II  »IMAX*  IOUlT.JLAST»N»Nl.N2»NBARl»NPiAR2t\lCYCL»MM\>HA»NN»NOUTP» 

2NPL3H,\OLHA«NPLUS*PI\'3»PINITtRADlU»SHPR»SIGMAtXBABl»XST0P,2MEXPt 

3T 

COMMON  APEA (200*2) .RURNE(2)  .CQSQXO) »CSQI200) »DTMI \ ( 3 ) »DTSU ( 200 ) , 
IE (700 *2) *FIVT(3) »EKI w ( 3 ) »F ( 200 *2 ) *HALFM( 700 ) .HALFR ( 3 ) ♦ INTER( 4 ) *  { 
2M20NE  (  3  )  *P  ( 200  *2  )  »PPLUS  (2  00).  *0(200*2)  »ROZER  ( 3  )  *  THETA  (  200)  »  ' 

3TVINS(2) *U( 200.2 ) *USO ( 200 ) »V ( 200 *2 ) »X ( 200 *2 ) 

1  CALL  LIRE 

1100  FORMAT ( 40X » 1 9H  INITIAL  CONDITIONS) 

1101  FORMAT ( 35Xt 15^  MASS  EXPLOSIVE.9X.F10. 2. 8H  GRAMMES) 

1102  F0RmAT(35X.17H  INITIAL  PRES5URE.7X .F10.4.4H  PSD 

1103  FORMAT  (3*jX,16H  MASS  PROJECTILE.8X.F10. 4. 8H  GRAMMES) 

1106  FORMAT ( 35X *22H  RADIUS  0.  THE  CHA.mBER*2X,F10.4»3H  CM) 

1107  FORMAT ( 35X» 22H  AREA  OF  THE  BARREL  .2X.F10.4.4H  C *2) 

1104  FORMAT (1H  ) 

110P  FORMAT (35X.20HPPOJ.  RELEASE  PRESS. .4X.F10.4.4HBARS) 

1109  FORMAT ( 35X.15HC0UNTERdRESSURE*9X*F10.4»4HT0RR)  | 

'/.'RITE  (6.1104 
U'R  ITE  ( 6  *1 104 
WRITE (6.1100 
MRITF(6.1104 
WRITE (6.110] 

WRITE ( 6.1104 
WRITE (6.1107 
WR ITE ( 6 • 1104 
WR ITE ( 6  » 1103 
WRITE (6.1104 
WR ITE ( 6 . 1 106 
WRITE (6.1104 
WR ITE ( 6 . 1 107 
'•'RITE  (6*1104 
VRITF ( 6  *1 10P 
WRITE (6.110A 
WRITE (6.1109 
CALL  DEPOT 

2  IF(BRA<E-1.) 

803  TF(PPLUS(JLAST-1)-SHPR)  801.802.802 

801  IMAX»2 
°RAKF“0. 

GO  TO  805  • 

802  IVAX-3 

ARA<F»1.  I 

805  CALL  GUTS  $ 

CALL  LFTOV 

IF(TOUIT-I)  2.6.6  I 

6  IF (DUO)  7.7.1 

7  STOP 

END  ’ 


) 

) 

) 

) 

ZMEXP 

)  PINIT 
) 

)  EM PRO 
) 

>  RADIU 
) 

)  AREA1 
) 

)  SHPR 
) 

)  PIN3  j 

803.802.807  l 


it; 


SIRFTC  READ  NOLIST 
SUBROUTINE  LIRE 

COMMON  A9EA1 *BRAKE»CSQyA*DtAMA*DTLA$*DUO»EINSU*EKSUM*EMPRO*ESUM» 

1 I » IMAX ♦ IOUI T * JLAST»N» Ml ,N2 »NBAR1 .NBAR2 »NCYCL.NMNHA»NN *NOUTP ♦ 
2NPL3H  *NPLHA  *NPLUS  *PJ!S3.jPJ  NIT  *RAQ1LI*SHPRaSIGMAaXBAR1  *XSTQP  tZMEXPt 
3T 

COMMON  AREA (200*2) *BURnE ( 2 ) .CQSOX ( 3 ) »CSQ { 200 ) *DTMlN ( 3 ) tDTSQ ( 200 ) » 
1E(200*2)#EINT(3) »EKIN(3) *F(200*2) »HALFM(200) »HALFR(3) » INTER (4) » 
2NZ0NE ( 3  J  *P ( 200 *2 ) *PPLUS ( 200 ) ,0 ( 200  * 2 ) *R02ER ( 3 ) , THETA ( 200 « * 

3TMINS <  2 ) *0(200*2 ) »USQ( 200 ) »V(200  *2 ) tX (200*2) 

1  FORMAT (313) 

2  FORMAT (8F10. 8) 

3  FORMAT (F10. 4) 

READ (5*1)  IMAX 

READ( 5 *2)  APEA1*Ra'OIU.XSTOP 

READ (5 *3)  EMPRO 

READ(5*3)  Z^EXP 

READ(5»3)  PINIT 

READ (5 *2)  SHPR 

READ (5*2)  PINS 

READ (5*1)  N ZONE ( 1 ) »NZOME ( 2 ) 

READ ( 5  *2  )  XRAR1 

READ (5.1)  NRAR1.NBAR2 

PEAD(5*1)  N1*N2 

READ (5 *2)  DUO 

RETURN 

END 


SIPFTC  SFTUP  NOLIST 
SUBROUTINE  nEPOT 

DIvFNSION  OUTBD (  2  )  .FZERO(  3  )  »V?FRO(  3  )  tOELXd  )  ♦GA^VA ( 3  ) 

COv'*ON  AREAl .BRAKE »CSOHA.DLAVA.DTLAS*DUO»ElNSU»EKSUV.£MPRO*ESUV» 

1  I » I HA  X . I OU I T  .JLAST  •  N .  ;•>  1  . \2  .  NBARl  »,\B AR2 » \CYCL  » NMNHA »  <v\  •  NOUTP  . 
2vpL3H  *MPLHA  *NnLUS  »PI  '-3  »P  I*’I  T  .RADIU  »SHPR  »SlG'‘iA  .XBAR1  .X5T0P  »Zv«tXP  * 
3T 

COvvON  AREA(200»2)  *euRNE(2)  »CQSOX(3)  *CSQ<200)  fDTMl.\<3>  »DTSa(200) » 
IE (200.2) .FI  NT (3) .EKl N ( 3 ) .F ( 200 »2 ) .HALFV ( 200 ) .HALFR ( 3 ) » INTER (4 ) . 
2vZONE ( 3 ) »P { 200 .2  ) * PPL US ( 200 ) .0(200.2) .ROZER ( 3 ) .THETA ( 200 )  . 

STAINS (2) .U(200.2) .USO( 200) ,V (200.2) »X( 200.2) 

RRAKE“C. 

COSOX ( 1 )«6. 

COSOX ( 2 ) >4* 

COSOX ( 3 ) *4. 

VZERO ( 1 )  «1» 

VZERO ( 2 ) *1 ♦ 

VZERO ( 3 ) “1* 

GAVYA(3)»1.4 
RURNE ( 1 ) »1 • 

°URNE ( 2 ) *1 • 

OUTRO ( 2 ) “RADIU 
R07F9 ( 1 ) »0.588 

ROZER(2)«1?./(14.5*300.*83.70)*PINIT 
ROZER (3)»28. 96/ <  750. 0*300. 0*83. 17) *PIN3 
FZERO< 1)«30400. 

FZFRO(2)-l. 225*1. E  05  *ROZER<2) 

EZ ERO ( 3 ) *P 1^3/ ( ( GANVA (3) *1.0) *7 50.0) 

VEXP«Z*EXP/ 0. 588*3. /2./3.U16 

CEXP«RA0IU**3-VEXP 

IF(CEXP)  1094.1099.1094 

1094  Yl«l, 

1095  Y2»Yl-( Y1**3-CEXP)/ (3.*V1*#2 ) 

IF(ARS( (Y1-Y2J/Y2)  -l.E-4)  1096. 1096 » 1097 
1097  Yl.vr 

GO  TO  1095 

1096  Cl)BRT«RA0  IU-Y2 
GO  TO  1030 

1099  CURRT  «0. 

1030  OUTRO ( 1 ) «CURRT 

1000  M*1 

1001  mv.MHA  «•" 

100?  'JPL3H 

1003  NPLUS  “1 

1004  VPLHA  “NPLUS 

1005  NN»3 

1006  MOUTP«0 

1007  NCYCL  “0 

1008  T-0. 

1009  OTLAS“0. 

1010  DTVIN(MN)«0. 

1011  INTER  (1)“1 
101?  DO  35  I«1.I»AX 

35  HALFR  (  ! ) “ROZER  (I)/2. 

N70VF ( 3 ) “NBAR1+NBAR2 
DO  40  I“1 *  I WAX 

40  INTER  (1+1) “INTER  ( I ) +NZONE  (I) 

JLAST* INTER (I  VAX) 


JLAS3* INTER { I VAX+1 > 

JLAS1*JLAST-I 

DNZOi*NZONE(l) 

DFLXUi*  (RAOtU**3-CRADlU-OUTBDtl) )**3)/DNZ0l 

VGAZ*  < KADIU-OIJTBD  ( 1 )  )  **3 

X(1»N)*0. 

U(1.NJ*0* 

00  55  1*1 tt MAX 

IF ( I .EQ.l.AND.NZONE ( 1 ) «EQ«  0 )  GO  TO  55 
JMIN-INTER  U)+l 
JMAX* INTER  C  I-^l  > 

00  51  J*JMIN»JVAX 
JMNHA  *J-1 
U( J*NMNHA)«0. 

0( JMNHA»N)*0. 

V( JMNHA  »N) *VZERO( I ) 

E ( JMNHA  #N)*EZERO(IJ 
IM  J.2)*G. 

Q( J^NHA.2)*0. 

V( JMNHAt2)*VZER0(I) 

E ( JMNHA » 2  J  *EZtRO ( I ) 

IF ( 1*2 )  101 .101.108 
C  SETUP  REGION  1  AND  2 

101  DL-J-1 

V0L*0ELX(1 >*3.14159*0.66666566 
HALFMt JMNHA)*HALFRl I >  *VOL 
CJ*RADIU**3-DELX  Cl) *D  L 
GO  TO  103 

102  PART*NZONE ( 2 ) /2 
M»ffZONE  ( 2 )  /2+NZO.NEJ 1 )  +1 
IF(J-M)  104.104.105 

104  DL*J-1-N20NE(15 
V0L*4./5.#VGAZ/PART*3. 14159*0.6666666 
HALFM ( JMNHA  >  *HALFR ( I ) *VOL 
CJ*VGAZ~4./5.*VGAZ*DL/PART 

GO  TO  103 

105  PART 1 *NZONE ( 2 ) /4 
MP*M+NZ0NE(2>/4 
IF(J-MP)  106.106.107 

106  OL*J-M 

VOL*3./20.*VGAZ/PART1*3. 14159*0. 66666666 
HALFM t JMNHA J *HALFR J I > *VftL 
CJ*1. 0/5. *VGAZ-3 .0/20. 0*VGAZ*DL/PART1 
GO  TO  103 

107  DL*J-MP 

VOL*l .0/20. 0*VGA7/PART1*3. 14159*0. 66666666 
HALFM ( JMNHA  > -HALFR  ( I ) *VOL 
CJ*i.0/20.0*VGAZ*( (PART1-0L) /PARTI) 

103  IF(CJ)  94.99.94 

94  Yl-1,0 

95  Y2*Y1“( Y1**3*CJ ) /( 3. *Y1**2 ) 

IFUBSnYl-Y2)/Y2)-l.E-4>  96.96.9  7 
Y1-Y2 
GO  TO  93 


97 


96  CUBRT  »Y2 

GO  TO  50 
99  CUBRT  *0» 

GO  TO  50 

C  SETUP  REGION  3 

10$  P(JVM1A»LI«PIN3/75Q* 

PI JMNHA»2>*PIN3/750« 

JDJ«J-JLAST 

IF(JDJ-NBARl)  109*109*110 

1 09  DELX1 *XBAR1/FL0AT(NBAR1) 

V0L*DELX1#AREA1 

HAUPV  CJ-^wwA  1  *vOL*HAjLF  B  LU 
DJ*JDJ 

CUBRT*-DJ»DELX1 
GO  TO  50 

110  IF (NBAR2 )  50*50*111 

111  XBAR2>l*2*XSTOP-(RADIU+XBARl) 
r>EUC2>X,BAii2 /FLOAT  (NBAR2) 

V0L*DELX2#AREA1 

HALFM ( JMNHA ) ■VOL*HALFR  f I ) 

0J-J0J-NBAR1 
CUBRT«-XBAR1-DJ#DELX2 
C  END  SETUP  REGION  3 

50  XU*N)«RADIlXUBRT 
51  X( J»2)-X(J*N) 

UUMAX+1*NMNHA)«0* 

E(JMAX*N)»EZERO(n 
55  V ( JMAX • N ) * VZERO tl) 

JLAS2-JLAS3-1 

HALFM { JLAS3 J iHALFMt  JLAS2J 

DO  10  J*1 » JLAS3 

J*J 

IF  (X  l  JfN'PLUS  J-RADIU  >400*5.;  0*500 
400  AREA! JtNPLUS  > ■ ( (RADIU  -X ( J.NPLUS  ))**2 >#?*#3.1416 ' 
GO  TO  10 

500  AREA(J*NP.LUS- J*AREA1 
10  AREA ( J*2 ) «AREA ( J  *NPLUS ) 

P( JLAS3.N)«PlN3/750. 

P ( JLAS3  *2 )*PlN3/750* 

DO  130  1*1  * IM4X 

IFd.EO.l.AMD.NZONEd  I.EQ.O)  GO  TO  130 
1*1 

IFU-2)  72*71*73 
72  CALL  EQST2 
GO  TO  130 
71  JMIN*INTER(I)+1 
JMAX*INTER(  Id  J 
JMIN!*JMIN-1 
jMAX1«JMAX-1 
F(JMAXl*N) *0*999 
F { JMAX1 *NPLU5 ) *0  *99n 
CALL  EQST1 
GO  TO  130 
73  CALL  EQ5T3 
130  CONTINUE 

DO  135  JMNHA*1 » JLaSS 
135  PPLUS  (JMNHA  )*P(JMNHA  *N) 

EINSU-0. 

EKSUM-O. 


FSUV-O. 

US0(1)«U(1»N)##2 
HO  165  I «1 • IMAX 

einti  n*o. 

FMN(  I  )*0. 

IF(I.EQ*l.ANDtNZONE(l) .EQ.O)  GO  TO  165 
J^IN* I NTFR  (I) 

JMAX«INTER  (I+1J-1 

no  160 

JPLHA  *J 

USO  ( J+l ) *U  (J+l  *N  )  **2 

FI  NT ( I ) “FI NT ( I  )+E (JPLHA  .NM/HALFR  (  I  )#HAl.FM  t  .IPLmA  ) 
160  EKIN( I >*(USQ( J)+USQ( J+l) )*HALFM( JPLHA  )  +EKIN ( I ) 
EKIN( I )*«5#EKIN( I ) 

EINSU  ■EINS'J  +EINT( I ) 

EKSUM*EKSUM+EKIN ( I ) 

165  CONTINUE 

FSUMbEIN.SU  +EKSUfc? 

CALL  SORT  I 

MPLUS  *2 

NPLHA  -NPLUS 

DTMIN(NPLHA  )*0*00006 

DTMIN ( NN )»Q* 00006 

TMTNS  (NPLHA  )«.00000008 

I0UIT-0 

X  ( JLAST  tNPLUS )  *X  <  JLAS  T.t.N ) 

RETURN 


SI«FTC  GKOS  \'OL  1ST 
SUBROUTINE  GUTS 

DIMENSION  DLAMD(200) »TSIGS(200) 

COMMON  APEAl»nRAKE*CSQMA*DLA‘<A*DTLAS*DU0*FI,\SU*EKSU.'.*EMPR0*E5UM* 

1 1 . 1 v.AX  ♦  I OU I T  » JLAST »N»N1  ♦  v2 » NBAR1  y\BAR2  ♦  MCYCL . NV.'iHA *.'JN  » NOUTP » 
2NPL3idji,‘^LHAjtJSl£l.USjPJ^.tP  I  M T  T  *R  AO  T  l.l  »  SHPR  i  S-IG  V.A  >  XBAR1  »X  STOP  »£X&XP» 

IT 

COTTON  AREA(2O0*2) »BURNE<2) »C0SQX(3) *CSO<200) »DTMlN<3) »DTSQ(200) ♦ 
IE  (  200  » 2  ) *EINT(3)  » FK I  'i (  3 )  »F  (  200  »2 )  *HALFV  (  200  )  »HALFR(  3  )  » INTER  (4)  » 
2NZ0NE ( 3 ) *P ( 200  *2 ) * PPL US (200  #0(200*21  .ROZER ( 3 > .THETA ( 200) » 
3TMINSI2)  ,U(?00*2  )  *'JSO  <  200 )  *V  { 200  *2  )  *X( 200*2  ) 

.iLAST«J^TE5(3) 

JLAS1-JLAST-1 
FRICOI.O 
OLAMA  -0* 

SIG’-'A  -0. 

T«T+DTVIN(NPLHA  ) 

VC.XCL  -NCYCL  +  1 
NOUTP-NOUTP+1 
DO  245  I  -1 »  I. VAX 

IFU.EO.l .AvO.NZGMEm.EO.O)  GO  TO  245 
CSQYA  -0. 

X  ( 1 »  N  )  «  0  ♦ 

X ( 1 »NPLUS) -0* 

U(1»NPLHA)*0. 

U ( 1 »NVNHA ) *0* 

I-I 

JK IN- INTER  <I)+1 
J^AX- INTER  <1+1 ) 

DO  196  JsJtflNt JMAX 
JPLHA  -J 
JVNHA  -J-l 

IF(JLAST-J)  806*807  #806 
807  IF(BRAKF-1.)  808*809.309 
80p  nUDT«0. 

GO  TO  195 

809  nuDT«FRICC*(PDLUS<J',VMA)-PPU'S<  JPLHA)  )*AREA1/ 

1  ( HALFM  ( J.XNHA )  +HALF  - ( JPLHA  )  +EMPRO ) 

GO  TO  195 

806  DUDT- ( PPLU5 (J'NHA  J-FPLUS (JPLHA ) ) *ARCA ( J»N ) /. 

1  <  HALFV  <  JVNHA )  +HALF’*'  <  JPLHA  )  ) 

195  I’U.NPLHA  ) *U <. J * .V-'NMA  )+DTVl!V(.\iN)*DUDT 
196  X  ( J.NPLUS  )  -X  ( J  *  *•! )  +DT>'  I N  ( NPLHA  )*U(J*NPLHA  ) 

JX|N* INTER! I )+l 
JVAX- INTER  C 1+1 ) 

C  AREAS  AND  VOLUMES  ARE  NOV.*  CALCULATED  FOLLOWING  THE  SCHEME 
C  OF  FLAGG 

00  230  J-JVIN.JMAX 
JPLHA-J 
•,'*.\,HA*J"1 

PRO-SORT (AREA1/3. 141 59) 

IF  CX { J»NPLUS)-(RADIU-PRO) )  400*500*500 
400  AREA< J»NPLUS)«<RADIU-X(J»\PLUS) )«*2*2. 0*3. 14159 
VOLUM-2#0*3 • 14159*< RAOIU~X ( J*NPI  US ) ) #*3/3*3 
GO  TO  17 

500  I F ( X ( J  *  NPLU5 ) “RADI U )  510*510*520 
510  AREA <J.  "PLUS >»3. 14159* <RADIU-X(J*NPLUS ) )**2+AR£Al 
IF(X( J-l »NPLUS)-( RADI U-PRO) )  512*512*513 
512  V0LUV«2*0*3. 14159* (RAD IU-X( J-l.NPLUS) ) **3/3.0 

l-3.l4l5o*(RAnTU-X(  J*NPLUS)  )*(  (RADIU-X1  J.N.PLUS)  )**2+3.0*PRO**21/6- 
GO  TO  17 


513  V0LUM«3.14159*( (RAOIU-XJ J-l.MPLUS)  >  #(  ( RADIU-X  ( J-l  .NPLUS  >  )**2+3.0* 

1PR0#*2)-(RADIU-X!  J.NPLUS)  )*(  (RADIU-X (J.NPLUS) ) **2+3.0#PRO*#2 )  )/6*0 
GO  TO  17 

520  AREA( J«NPLUS)*AREA1 

I F ( X ( J-l. NPLUS J-RADIU)  523.522.522 
5??  V0LUV«AREA1*(XU.NPLUS>-X(J-1.NPLUS)  ) 

GO  TO  17 

523  IF (X ( J-l.NPLUS)-(RADIU-PRO)  )  524*524.525 

524  V0LUV«AREA1*(X (J .NPLUS J-RADIU) +2. 0*3*14159/3. 0* ( RADIU-X ( J-l .NPLUS ) 
1)**3 

GO  TO  17 

525  V0LUM«AREAHMX(J.NPLU5)-RADIU)+3. 14159/6. 0# ( RAD IU-X ( J-l.NPl  US) ) 

1*( (RADIU-Xl J-l .NPLUS) )*#2+3.0*PR0*#2 ) 

17  CONTINUE 

V( JMMHA  .NPLUS  )-VOLUM  /  HALFM ( JMNHA  )  *hALFR(I> 

IF(U( J.nplHA  )-U ( J-l.NPLHA  ))  205.225.225 
205  Q l JVNHA  .MPLHA  )»COSOX  (I)  *(U( J.NPLHA  )-U ( J-l .NPLHA  )) 

1**2/(V( JVNHA  .NPLUS  )+V(J«NHA  »N))  *HALER(D 
GO  TO  230 

225  OUMNHA  .NPLHA  )«0* 

230  CONTINUE 

IFU-2)  83.82.81 
83  CALL  EOST2 
GO  TO  53 
82  CALL  EOST1 
GO  TO  53 
81  CALL  EOST3 

JMIN-INTER( I )+l 
JMAX- INTER ( 1+1 ) 

53  00  240  J-JMIN.JMAX 
JMNHA  *J-1 

TSIGS  (JMNHA  )-C5Q(JMNHA  )/(X(J*NPLUS  )-X( J-l  .NPLUS  )  )**2 
PPLUS  (JMNHA  )-P( JMNHA  .NPLUS  )+0( JMNHA  .VPLHA  ) 

DLAMD ( JMNHA  ) -COSQX (I )/2«*(V( JMNHA *N)-V( JMNHA. NPLUS ) ) / 

1(VC JMNHA  »N  )+V ( JMNHA  .NPLUS  )) 

IF C DLAMD ( JMNHA ) -DLAMA )  600.600.601 

600  DLAMA-DLA.MA 
GO  TO  620 

601  DLAMA-DLAMD( JMNHA) 

620  IF(TSIGS(JMNHA)-SIGVA)  700*700.701 

700  SIGMA-SIGMA 
GO  TO  720 

701  SI GMA-TSIGS( JMNHA) 

720  IF(CSQMA-CSO( JMNHA) )  800.800.801 

800  CSQMA-C5Q ( JMNHA ) 

GO  TO  240 

801  CSOMA-CSOMA 

240  DTSQ(JMNHA  )«.1111111/TSIGS  (JMNHA  ) 

245  CONTINUE 

!F(XSTOP-X( JLAST .NPLUS) )  340.340.246 
340  I QUIT-1 
CALL  SORT I 
246  RETURN 
FND 


waggig«ggB3g5gSS55igl 


SleFTC  RETCUR  NOLlST 
MmwouiIME  LFTOV 

COWSN  A5EA3  f  r*?.AKE  »CSQVA  *DLA''A  .DTLAS  »DUO  *EI  NSU  .EKSUV.EVPkO  » ti»UM» 
1 1 »  I-’AX  »  IOUI T» JLAST»M*  Ml *N2  .f'PAftl  ♦NBAR2»NCYCL»N*/ihHA.NN  *NUU  I K  * 
2>PL*H,\>DLMA»Nr>LUS.P!r.3»PI.'!!TtRAD!U»SHPR»SIGMA*XRAPi  •XSUltuxMJtAK. 
3T 

CO'-  •■«,*!  AREA  1 200»  2  )  *P'JPNE  (  2  )  *  COSOX  (  3 )  *CSO  ( 200 )  » DTK IN (3  )  *0  1 6U(*00 ) 
IF ( 200  *2  )  *EI  \'T  (  3 )  »EKlf>  (  3 )  »F  ( 200  .2 )  »HALFM (  200  )  .HALFR (  3 )  »INIfcK(4)  » 
2f!?.0vEO)  »n>  (?00»2  )  »PPLUS  (200)  »Q  ( 200  *2 )  »R02ER  ( 3  )  #THETA(  200 )  * 
3T;-'I\S(2)  »U  (200*2  )  »USO(?OQ)  » V  ( 200  *2  )  *X  ( 200  *2  ) 

JLAST* I NTER ( 3 ) 

DTLA'.  *DTVI  -V  ( NPLHA  ) 

SIGVJ  -l./SIGNA 

I F  ( S I GM I  -9.  *TMINS  (NPLHA  ))  270*270*250 
250  IF(DLAVA  -.08  )  255*270*270 
255  IF (SIG^I  -16.  # TV I NS  (NPLHA  ))  265.260*260 

260  If  IDLAKA  -.05  )  270*270.266 
265  TN'INS  d’PL3H  ) -TWINS  (NPLHA  ) 

GO  TO  285 

270  IF(DLAWA-.0O5)  275.275.280 
275  TWINS  (NPL3H  ) «S IGMI  /4. 

GO  TO  285 

280  IF(SIGMI/16.  -.005184#TMlNS(NPl  HA) /DLAWA#42J  400. 400. 401 
400  TVINS ( N°L3H ) "SIGN*  1/16. 

GO  TO  285 

401  CM I MS ( NPL3H ) «. 005184*TKI NS (NPLHA ) /DLAMA##2 
285  I F  ( SORT  (TWIN’S  ( N°L3H)  )  -1»4*DTV.I  N  ( NPLHA ) )  500.500.501 

500  DT.VI  N  ( NPL3H )  -SORT  ( T-V.I  NS( NPL3H )  ) 

GO  TO  290 

501  DTVIN  ( NPL3H ) -1 »4#DTKI M (NPLHA ) 

250  DTMIN(NN)«(DTMIN(\PL3H  3+DTK IN (NPLHA  )>/2. 

NP-NPLUS 

NPLUS  -N  1 

NPLHA  -NPLUS 

M-MP 

NMNHA  »N 
MPL3H  -N 
FINSU  -0. 

ESUM-0. 

EKSUM-0. 

USO(l )-U(l*N)**2 
DO  300  I-l.IMAX 
EINTd  )-0. 

EKIN(I)«0. 

IF i I .EQ.l.AND.NZONE ( 1 ) .EO.O )  GO  TO  300 
JMIN-INTER  (I) 

JMAX- INTER  d+D-l 
DO  295  J-JWIN.jmaX 
JPLHA  -J 

USQ(J+1)*U( J+1*N)##2 

EINTd  ) '•El  NT  (I  )+E(JPLHA  »N)  ♦HALF?'!  (JPLHA  )  /HALCK(l) 

295  EK I N C I  )-(USO( J)+U5Q( J+l)  )*HALFM( JPLHA  )  ♦EKINU) 

EKIN( I )-.5#EKlN( I J 
EINSU  -EINSU  +EINTC I ) 

EKSUM-EKSUM+EKIN ( I ) 

300  CONTINUE 

EKSUM-EKSUMf.5#EMPRC  *U( JLAST.N)*#2 
ESUM-EINSU  +EKSUW 


IF(NCYCL-Nl)  120*324*324 
324  IF (NOUTP-N2 )  120*327*327 
327  N0UTP*0 
CALL  SORT! 

120  RETURN 
END 


P( J+l/2 *N 
DM( J+l/2) ) 
BARS 

GRAMS) 


P( J+1/2*N 
TEMPR  ) 
BARS 

DECK) 


S1PFTC  OUTPUT  .SOL  I  ST 
SUBROUTINE  SORT  I 
nr-'F.NSlO'J  2VASS (  200  ) 

CO  WON  AREA 1 »BRAKE ♦  CSOMA  *  DLAMA  »DTLAS  » DUO  »  E I NSU » EKSU.M  *  EMPRO  *E5UM  • 

1 1  j  IMA X  U OUT  T  Wl  A «  M  »  M |  ffoP  «  MRAR 1  >NftABg.^irVf  1  . .VM.VH  A  . MM  .  Mftl  l,TB . 

2NPL3H»NPLHA»NPLUS»PIN3»PINIT»RADIU»SHPR»SIG1’viA»XRARl»XSTOP»2MExPt 

IT 

COMMON  AREA(200»2) »BURNE(2) *C«SQX(3)  »CS0(200)  *DTV  IN  ( 3  )  »DTSQ  ( 200 ) » 
IE  ( 200  *2 )  »E  I  NT  {  3 )  *EKIN(3)  »F  ( 200  »2 )  »HALFM{  200 )  *HALFi<  ( 3 )  »INTER  (  4 )  » 
2NZONE (3) *P ( 200  *2  ) »PPLU5 ( 200 ) *0(200*2) *R02ER ( 3 ) * THETA ( 200 ) » 

L2J  a!1  L2G0 «2J  «LLSQ(.?(IQJ  »V U&U2 i,U2QiU2  ) 

IMAX«3 

JLAST* INTER ( I MAX  5 
JLAS1-JLAST-1 
JLAS3*tNTER  ( IMAX+1 ) 

JLAS2*Jl.AS3-l 
421  FORMAT tlHOJ 

1  FORMAT (1H1) 

2  FORMAT  (3X»115HJ  X(J»N)  U(J»N-l/2)  V(J+1/2*N)  P(J+1/2*I 

1)  Q(  J+l/2  *N-l/2 )  E( J+l/2 *N )  AREA ( J*N )  DTSO ( 1/2 »l/2 )  DM(J+l/2)) 

3  FORMAT (3X.111H  CM  CM/MJLLISEC  CC/CCO  BARS 

1  BARS  BARS-CC/G  MILLISECSO  GRAMS) 

A  EQflMAT  (_IiL*  6E.1 1.  5  »  2£_LL*  3  *  FlA*Al 
6  FORMAT ( 15  »  3E15.5) 

f  FORMAT 1 3X » 1 16HJ  X(J*N)  U(J*Nrl/2)  V(J+1/2*N)  PU+1/2* 

1)  0(J+l/2»N-l/2)  E(J+1/2»N)  AREA(J*N)  DTSO (1/2* 172 ) •  TEMPR  ) 
P  FORMAT (3XilllH  CM  CM/MILLISEC  CC/CCO  BARS 

1  BARS  BARS-CC/G  .MILLISECSO  -DECK) 

12  FORMAT  (47H  CYCLE  T  DT  TOTAL  Ei 

ILIMI  «IMAX 
WRITE (6*1 ) 

WRITE (6*12 ) 

WRITE (6*6)  NCYCL  «T »DTLAS  *E5UM 
WRITE ( 6  *421 ) 

IF (NCYCLJL  53 t.SL3>S4 
53  DO  57  JPLHA*1 * JLAS2 
57  ZMASS ( JPLHA ) *2  *#HALFM ( JPLHA ) 

WRITE (6 *2) 

WRITE(6*3) 

WRITE (6*421 ) 

DO  604_L"ltJLAS2 
IF(L-INTER<3) )  601*603*602 
601  IF (L-INTER ( 2 ) )  602*603*602 

602  WRITE (6*4)  L*X(L»N) »U(L»NMNHA) *V(L»N) *PPLUS(L) »Q(L*NMNHA) • 

1E(L»N) »AREA(L*N) *DTSQ(L) *ZMASS(L) 

GO  TO  604 

603  IF{L*FQ*LNTER(2)  .AND.NZ0_NEUJ.*£a«0)  GO  TO  6JQ2 
WRITE (6*421 ) 

WRITE (6*4)  L*XIL*N)*U(L*NMNHA  ) »V.(L*N) .PPLUS  (L)»Q(L*NMNHA  )* 
1E(L*N) »AREA(L*N) »DTSO(L) *ZMAS5(L) 

604  CONTINUE 

4 1  FORMAT ( 14*  2E13.5*74X*E13.5) 

WRITE 1 6  *41  )_JLAST *X ( JLAfiTtN l.u( JLASTtNMNHA) .EMPRO 
GO  TO  71  9 

64  IF (BRAKE-1*  >  80*81*81 
60  L2«JLAS1 
GO  TO  65 


U(J*Nrl/2)  V ( J+l/2 *N) 

AREA ( J*N )  DTSO (1/2 *172)-  1 

CM/MILLISEC  CC/CCO 

.MILLISECSO  t 
DT  TOTAL 


81  JL1* JLAST+1 

00  82  JaJLl* JLAS3 
JMNHA* J-l 
APIN«2.*PIN3/750# 

IF ( PPLUS ( JMNHA) -API N)  83*83*82 
83  L2*J 
GO  TO  65 

82  CONTINUE 
65  WRITE (6*7) 

WRITE (6*8) 

WRITE ( 6  *421 ) 

DO  704  L«1 *L2 
ILIMI-IMAX 
DO  702  lL«le.IUMI 
IF(L-INTERdL) )  702*703*702 
702  CONTINUE 

705  WRITE(6*4)  L*X(L*N)*U(L*NMNHA) *V(L»N) >PPLUS(L) *Q(L»NMNHA) * 
1E(L*N)  *AREA(L*NJ  *DTSQ.(L)  *THETA_LLJ 
GO  TO  704 

703  IF  {.IJL.A.ECU2*  AND»NZQNE(.l.)  •  EGh.OJ  GO  TO  705 
WRITE (6*421) 

WRITE (6*4)  LiJULiN) »U(.L*NMNHA)  »V(L*N)  •  &RLUS(L)..*Q(L*NMNHA>  • 
IE (L*N ) *AR£A(L*N) *DTSQ(L) *THETA(L) 

704  CONTINUE 

WRITE (6 *4)  JLAST »X{ JLAST *N ) »U( JLAST .NMNHA) 

71  CONTINUE 
RETURN 
END 


SIBFTC  ETAT1  NOLIST 
SUBROUTINE  EOST1 

CON  VON  AREA1 » DRAKE  *CSQMA  *DLAMA  »DTLAS  *DUO  »EI  NSU  »EKSUM  »EMPRO  »E5UM  ♦ 
II ,IMAX»lOUIT»JLAST*N»Nl*N2»NBARi»NBAR2»NCYCL»NMNHA»NN*NOUTP» 

2 NPL3H  ♦  NPLHA  . NPLUS  »P I N 3  »P Lit  T  »R ADJUaSMPR  . SlfiMA  *XRAR  1  •  XSIX1P  »ZM£JtP  • 
3T 

COMMON  AREA(200*2) »BURNE ( 2 ) »CUSOX ( 3 ) »CSQ < 200 ) iDTMIN ( 3 ) »DTSG ( 200 ) , 
IE  ( 200,2  )  »EINT { 3 )  *EKI<M(  3 )  *F ( 200 ,2 )  »HALFM(  200 I  »HALFR( 3 ) *  INTER (4 ) » 
2N70.NE  (3)»P{200»2)  »PPLUS( 200 ) »Q ( 200 ,2  )  »R02ER  1 3 1 ,ThETA( 200)  » 
3TMINS12) ,U(200,2 ) ,USO ( 200 ) ,V (200,2) *X ( 2uO,2 ) 

I  *2 

JY IN- INTER (I >  +  l 
JVAX*INTER( 1+1) 

JMIN1-JKIN-1 
J.YAXl-JMAX-1 
IFfBURNE  (I))  1,10,1 

1  FV1N-2, 

IF(NCYCL)  500,500,501 

500  JMAX2-JVAX1-1 
GO  TO  502 

501  JWAX2-JYAX1 

507  DO  66  JMNHA*JVIN1 » J'MX2 

F  ( Jf‘NHA  ,MPLUS  )»(1a«V(.JMNHA  .NELLIS  ))/(l#-‘  54) 

I F ( F ( JVNHA ,NPLUS ) “•00l'2,3,3 

2  F(J^NHA  .NPLUS  )«0. 

3  I F ( F ( JMNHA  ,NPLUS  I-l.)  4,6*5 

4  IF(F(J.v.NHA  „NPLUS  )-F(JMNHA  ,N))  5*6,6 

5  F ( JVNHA  .NPLUS  )«1. 

6  IF(F( JMNHA, NPLUSI-FMlN)  300*300,301 
“300  F V| I N- F  ( JMNHA , NPLUS I 
GO  TO  66 
301  FVIN-FMR* 

66  CONTINUE 

IF(FMI,N-1.)  10,7*7 

7  DO  8  JMNHA  *J--IN1,JMAX1 

8  FtJVNHA  , N ) * 1 • 

F ( JMAX1 *N)-1. 

F( J“AX1*MPLUS)»1# 

PURNE  (II-G. 

10  00  200  JMNHA  «JMIN1»JMAX1 

El  *  ( E  ( JV’JHA  *•')-(  (PCJMNHA  ».\)+Q(  JM-MHA  *NELbA  >># 

1  ( V  (  JmvhA  ,NPLUS  )-V(JMNHA  »NM  )  HUOOOl/ROZER  (II 
VE*V(JVNHA  *NPLUS  I/ROZER  (I) 

IF(F(J-VVHA  ,NPLUS  I-l. 0)  90,99*99 
9P  01-10. 

GO  TO  101 

99  Pl-P(  JMNHA  ,N  I4.000JL 
101  PlVE-°l#VF 

I F (P1VE-1 ,0465  I  11*11*12 

11  ALPHA-0.0 
DALFD  -0.0 
DALFV  -0*0 
GO  TO  15 

1?  IF(Jlvr-3.4«8)  13,13*14 
13  ALPHA«6AO;),0#P1VE-9000.0 
DALFD  -8600. 0#VE 
DALFV  -8600. 0#P1 


GO  TO  15 

14  ALPHA-21.0E03 
DALFD  *0.0 
DALFV  *0.0 

15  r>lvE2*PlVC**2 
P'J  VF4*PlVE2**2 
P1VE5«P1VF.4*P1VE 
PPAC1*1140.0+P1VE4 

r>PAC?«.10lF-03*ALOG(  ARS1P1/1.013E-03)  ) -. 232 5E-03 
EF*6.57#PM>974.0*P1VF2/RRAC1-ALPHA*BRAC2“E1 
EPRl'-'l  =6.57*VC+1948.0*VE**2*P1/BKAC1-3B96.0*VE*P1VE5/BRAC1*#2 
1-QALEO  #P2AC2-ALPUA*.101£-03/P1 
PNEV*Pl-EF/KPRI.v 
IF ( ABS (  PNEW-P1 )  -  .0001  )  17.17.16 

<n«D.\j£w 

GO  TO  101 

17  E  ( J'-'MHA  .NPLUS  )*El-(  . 5* ( P1*10C00. *F (  JYNHA  .KPLUS  )-P(jyNHA  »N))* 
llYU'^HA  .M.PLLLS  J“Y  1  JMtJfclA  «M  )  U.OOOl/RQZJiil  l!J 

171  PlVE*°l#VE 

PlVE2*PlVE*#2 
PlVE4«PlVE2**2 
PlVE5“PlVE4*PlVE 
IF(P1VE>1.0465)  18.18.19 

18  ALPHA*0.0 
DALFD  *0.0 
DALFV  *0.0 
GO  TO  2? 

19  IF (P1VE-3.488 )  20,20,21 

20  ALPHA*  8600  «*P1VE*"9000.0 
DALFD  *86fiajL0*V.E 
DALFV  *8600. 0*P1 

-G6  TO  22 

21  ALPHA*21.0E03 
DALFD  *0.0 
DALFV  *0.0 

22  a8ACl«U40«Q*PlVEA 

RRAC2«.101E-03*ALOG(ARS(P1/1.013E-03) I-.232&E-03 
£F»6,57*PlVE+974.0*PlVE2/BRACl-ALPHA*BRAC2-tC  JX.NHA  .NPLUS  ) 

FPRI'a  *6. 57*VE+I948.0*VE**2*Pl/BRAC 1-3896. *VE#P1VE5/BRAC1**2 
1 “DALFD  *PRAC2-ALPHA*. 101E-03/P1 
PMEW-P1-FF/EPPIM 
IF (.AflSJ  P'JEW-Pl  )r_«0001 )  24.2A..23 

23  Pl-PNEW 
GO  TO  171 

24  .E(J*NHA  .NPLUS  )aE(JVNHA  .NPLUS  )#1.0E04*ROZER  (D 
p ( Jv.NHA  .NPLUS  >*P1*1.0E04*F(J.v1i\;hA  .NPLUS  » 

PEDP*EPRI>/ 

0F0V»6  «  57»P1  +  194B«0»P 1VF»P 1 /BRACl-3fi96.0»PlVF5»Pl /8RAC1»»2 
1-DALFV  *BRAC2 

IFIFUMNHA  .NPLUS  >>202,201,202 

201  csoijMnha  )*1.0 

GO  TO  200 

202  CSO ( JMNHA  )«IV<JMNHA  .NPLUS  J/ROZER  (I )) **2« ( Pl+DEDV ) /DEDP*i.E04 

200  THETA ( JVNHA  ) *P <  JMNHA .A'PJLLIS ) «V ( JMNHA iNPLUAJ /  (A5 .  *R02F8LLU  * 12 # 

RETURN 

END 


SIRFTC  ETAT2  NOLIST 
SUBROUTINE  EOST2 

Convey.'  AREA  1  .BRAKE  »CSQyA*DLA.MA»DTLAS»DUO»  El  NSU*£KSUM»EMPftO»ESU'vi* 
II  »I''‘AX.laUIT*JLAST.N»Nl*N2*NBARl..\BAR2»NCYCL.NM,;HA»NN*NOUTP» 

2NP1  3H  #Al£L«A  «NPl  1LS  *R-1  N  3  *P  I  A!  II  >  S  AD  IU  »  SHPR  >  S  Ifi.ViA  »  XBAR I  »XST£P  .Z-afcXP. 
3T 

COV'^ON  AREA  (200*2  )  »BUR’NE<2)  *CQSQX(3)  *C$Q(200)  »DTMtN(3)  »DTSU(200) 
IF (200. 2)  ,EINT(3)  .EKINO)  *F ( 200.2 ) *HALFM( 200 ) »HALFK ( 3 )  ♦  INTER  14 )  » 
2NZ0NE ( 3 ) .PC  200 *2 ) .PPLUSC 200 ) *Q C 200 *2 ) *R02ER l 3 ) »THETA C 200 ) • 
3TMIVSC2) »U (200*2 ) *USO<  200) * / (200 *2 ) *X C 200 *2 ) 

1*1 

j\’l v* I \TFP  (  I  )+l 
JMAX* INTER  C 1  +  1  ) 

JVI\1«JVT\-1 
Jv AXl * JVAX—1 
IF  (BUR'-'E  (I))  1*10.1 

1  £yjt;«2-. 

no  66 JM,NHA  *JVIN1*J“AX1 

F  ( JV.NHA  .NPLUS  ) « <  1  .0-V  (  JKNHA  »N»LUS )  )  /  ( 1.0-0 .78  72) 

i  f  c  f  c  j'-'nma  » ''Plus  )  -. oo i )  2*3*3 

2  FCJMNHA  .NPLUS  )»0. 

3  IF(F(JvnwA  *MPLUS  )-l.)  4*6*5 

4  I F  ( F  C  J'-IMHA  »f<PUlS  )-F(JMNRA  »M)  )  5.6*6 

5  FCJvnhA  .'-PLUS  )»i. 

6  I F ( F ( JVNHA  » \ PLUS ) — Fy I N )  300*300*301 

300  FVI.\*F<  JWMHA»'!PLUS) 

GO  TO  66 

301  FMIN«FVI'' 

66  CONTINUE 

IF  CF'-'lN— 1.  )  10*7.7 

7  DO  8  JVNHA  *J.'  IN1.JVAX1 

8  FCJVNHA  *M)*1, 

PURNE  (I)*0. 

10  DO  200  J"NMA  «J' IN1.JMAX1 

F1»(E(J«'!HA  *C>  )-(P ( JiNHA  »N)-tO(OMAlHA  iNBLRA  )  U<VUA£liA  .ttPUlS  )- 
1 V ( J'-'NhA  *••■)  )  )*l.OE-06/R0ZER  (I) 

RHC*ROZER  ( I )/V( J^NHA  tNPLUS  ) 

RH0?«RH0**2 
EXPON*EXP (— 6./RH0) 

20  A«  .OC?164#vhO?*RHO?+;?.0755*EXPON 
P«*35*Shq 

OAORH  *  •0C°65S#RHO#'?HO2  +  12.45 30*EXPO\/RHO2 
DBDRH  *.35 
Pl*(A+B*Fl  )*1.0E06 

108  E ( J'NNHA  .NPLUS  )  *E1*1  •0E86*RGZS»  C  I  ) -.5* ( Pl-P ( J'-NHA  *N))# 

1  ( V  ( JVNHA  *  ’’PLUS  )-V(Jv.NHA  »v)J 

p ( JM'*HA  .'-PLUS  )«(A»t'»E.D>-lHA  *>iPLUS  ) *XaC.fc'-Gi> /liwZER  11)  )  *1.0£U6 
1*F(JV.\HA  .NPLUS  ) 

IF(F(Jv,vHA  *W>U/S  )  )202*?01*202 

201  CS0(J''Nma  )  *1.0 
GO  TO.  200 

202  CSO ( J*’f  HA  )«(H*P(  J-.’\'HA  .NPLUS  ) #1  .OF— 06/RM02  'DAD'-JH  ♦9HI.'Ph*‘ 
1E(J'--*'HA  .NPLUS  >*1.0ErO6/ROZFR  C I )  ) *  1  .££.06 

200  THETA  (  JvmmA  )*P(  J  *'.NHA*,vPLUS)*V(  JV\'HA*\PLUS)  /Ct!3.l44#^oZf:R(  I )  ) 
RETURN 
E*D 


SlBFTC  ETAT3  NOLIST 
SUBROUTINE  EQST3 
REAL  MU  »VU0  tf'U2 

COMMON  AREA1.9RAKE.CS0MA.DLAMA.DTLAS.UU0. El  NSU.ElCSUM.EMPRO.ESUM. 
II  » IMAX. IQUIT  •  JLAST *N»  N1.N2.NRAR1  .NBAR2 .NCYCL.NMNHA.NN .NOUTP. 
2NPL3H*  NPLHAJ  NPLUS  •  P  T  N  ^  » P I N  T  T  «  RAD.I  UaShPRaSIS^A  OtBARl OCSIXJP  »ZUiXP » 
3T 

COMMON  AREA  1 20C  *2 ) *BURNE ( 2 ) tCOSQX ( 3 ) »CSO ( 200 ) .DTMlN ( 3 ) .DTSQ { 200 ) • 
1F(?00»2) »EIVT(3) tEKIN(3) »F(200»2) »HALFM(200) »HALFR(3) »INTER(4) . 
2NZONE ( 3 ) »P ( 200  *2 ) .PPLUS ( 200 ) *0 ( 200  *2 ) »R02ER ( 3 ) »  THETA ( 200 ) » 

STMIMSt  2) »U(200»2 ) *USQ ( 200 ) »V (200*2 ) »X( 200.2 ) 

1-3 

JMlN*INTER(n+i 
JMAX*INTER( 1+1) 

JMIN1*JVIN-1 
JVAXiijMAX-1 
DO  >200*  JMNHA*JMIN1#JVAX1 

£!>(£( JMNHA »N>-<  ( P  t  JMNHA  >N  >  40  (  JMNHA  .NPL.HA  )  )  *  i  SJ  t  JMNHA  ...PL  US i 

i-»V(  JMNHA  »N )  ) ) ) /( 1960.0*ROZER  (I)  ) 

P1?P(JMNHA.N»/1.01375 

rtETAiROZER  ( 3 )  /  ( V  ( JMNHA  *\'PLUS )  *1  •  293F-03 ) 

PSI*ALOG(HETA) 

ZIfPI/HETA 

NE«1 

1  fpliO/Zl 

A»(25  #8949.684#Y+3 .0 )  /  ( 4.  778  97369*Y+1 .0 ) 
^DADY6ll*5579473/(4*77e97369*Y-fl.0)**2 
*  J:B*®6iiiOfY.#(l*0-fYJ/(3000.;)#Y*#2‘t-1.0) 

DBDY*ft61.0#(1.0-2.0#Y)/(3000*0*Y*#2+l.C> 
i-\0.5166E+07#Y##2*  1 1  .  0-Y )  /  ( 3000  .0#Y##2+1 . 0  l  ##? 

. C?Z536. 0#Y# ( 1 . 0-Y  > / ( 0 . 9E+05*Y**2+ 1.0) 
y:0CC^F2336«O#(  1.0<-2. 0#Y )/ <Q.9E+05#Y  **2+1.0) 
'iM*^56i48E+09fY##2*(l#0-Yj/(0.9E+05#Y##2+1.0)*#2 
^D«0i;4iE+654 ( 1# OrY )*Y/  ( 0. 12E+08*Y##2+1  *0 ) 
bDDY*0i4lE+05* (1 • 0-2 • 0#Y ) / ( 0 • 1 2E+0B#Y*#2+1 • 0 ) 

Vr0.984E+10#Yf»2» 1 1,.0-Y )  /  CD .12g+0fl»v»»2+l.n ) *#2 
MU0*i.0+A+9+O*D 
:DMUOYibADY+DBDY+DCDY+DDDY 
MU2)S.{  6002 . 0*Y+4, 0 )/( 1 000 .0#Y4l .  0  ) 

DM.U2Y  *2002.0 /( 1000  *0#Y+1 .0 )  ##2 
MUiMUO+0 .  09#  (  .VU0-VU2 )  #PS I 
DMU Y«,U  •  0+0  •  0,9*PS  1 1  «D MUOY-0  cQ9*PS I  *DiiU2Y 
bMU4-DMUY*Y-**2 
EF*0. 20#Z1* < MU— 1 ,0 ) -E 1 
'  EPR I M*0 . 2  0# ( MlJ-1  •  0+Z  1  *D."U ) 

ZNEW*  Zl-EF/EPRIM 
IFIARS ( ZNEW-21 ) -.0000 1 )  3.3.2 
2  Z1.PZNEW 
GO  TO  1 

3  I F:I NE-1 )  4.4.5 

4  v.R  1  ■.  1„.  0 1 3 7 5 #  H  E  T  A#  Z 1 

:El*El- ( 0.5#  < ( Pl-P ( JMNHA  .  N ) )  #  ( V  C  JMNHA  .NPLUS )  -V  { JMNHA  *N  ))> )  / 1 1960.0 
lfROZER (I)) 

NE*2 
60  TO-  1 

5  E:(  JMNHA  .NPLUS)  *1960  .0#R0ZER(  I  )#E1 

TJMNHA .NPLUS) *1.01 37 5#HETA#Z1 
3*485 «0/( 1000«C+21##2 1+3860.0  /( 7500.0+16. 5#Zi> 

THETA ( JMNHA ) *2  73 . 2#Z 1 *S 

200  CSO(  JMNHA  I  «1  *40#P  (JMNHA  .(jPLUS )  *V  ( JMNHA  «2iPLUS}  /RQZES  ( 1  i 
RETURN 
END 


SUBROUTINE  LIRE 


This  subroutine  reads  the  necessary  data  to  start  the  program.  There  are  12 
input  data  cards. 

t 

-  1st  card  IMAX:  has  to  be  taken  as  3  in  the  present  version  of 

the  program;  no  other  choice  is  possible 

2 

-  2nd  card  AREA1:  area  of  the  barrel  in  cm  . 

RADIU:  radius  of  the  barrel  in  cm. 

XSTOP:  length  at  which  the  program  is  to  be  stopped. 

-  3rd  card  EMERO:  mass  of  the  projectile  in  g. 

-  4th  card  ZMEXP:  mass  of  explosive  in  g. 

-  5th  card  PI  NIT:  initial  pressure  of  the  gas  mixture  in  PS  LA 

-  6th  card  SHPR  :  pressure  at  which  the  diaphragm  bursts  in  bars. 

-  7th  card  PIN3  :  initial  pressure  in  the  barrel  in  Torr. 

-  8th  card  NZONE(l):  number  of  zones  in  region  I  (I  =  1  and  2  ). 

-  9th  card  XBARI:  length  of  the  barrel  divided  in  NBARI  zones. 

-10th  card  NBARI:  number  of  zones  in  the  first  part  of  the  barrel. 

NBAR2:  number  of  zones  in  the  remaining  part  of  the  barrel. 

-11th  card  Nl:  cycle  number  of  the  first  output. 

N2:  cycle  interval  between  two  outputs. 

-12th  card  DUO:  if  DUO  <  0,  program  is  to  exit  after  con^letion, 

if  DUO  >  0  another  run  is  to  be  done 


Set  of  sample  data. 


CALL  LFTOV 


SUBROUTINE  DEPOT 


The  initial  conditions  for  the  problem  are  set  up  in  this  subroutine 


$ 


SUBROUTINE  GUTS 


This  subroutine  calculates  new  variables  for  each  time  step 


NCYCL 


SUBROUTINE  LFTOV 


New  time  steps  are  determined  and  total  energy  is  calculated 


SUBROUTINE  SQRTI 


This  subroutine  provides  printed  output  when  print-out  cycle  is  reached. 


Ji 


SUBROUTINE  EQST1 

Equation  of  state  for  gas  mixture  2H^+  0,,(see  Sec.  2.3.1) 


SUBROUTINE  EQST2 


Equation  of  state  for  explosive  (see  Sec. 2. 3. 2) 


gg"' . . 


SUBROUTINE  EQgT3 


